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PREFACE 


This  study  of  the  motion  of  earth  satellite  vehicles  in  trans- 
lunar  space  was  undertaken  at  the  request  of  the  Directorate  of  Space , 
DCS/Research  and  Development.  While  it  is  sometimes  possible  to  es¬ 
tablish  a  satellite  about  the  earth  in  a  translunar  orbit  that  will 
remain  reasonably  circular  for  as  long  as  five  years,  it  is  desirable 
to  have  an  orbital  correction  capability  tc  reduce  orbital  variations 
caused  by  the  sun  and  the  moon.  In  this  report,  a  method  is  described 
for  determining  orbital  injection  conditions  that  minimize  the  devia¬ 
tions  of  the  orbit  from  circularity  and  a  method  of  orbital  correction 
is  presented  that  can  be  used  to  reduce  variations  from  circularity 
if  they  develop.  The  results  obtained  and  the  methodology  developed 
herein  should  be  of  use  in  predicting  and  controlling  the  orbital  mo¬ 
tion  of  the  satellite  vehicles  in  regions  where  the  gravitational 
perturbation  due  to  the  moon  and  sun  are  significant  in  comparison 
with  the  earth* 8  gravitational  attraction. 
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SUMMARY 


This  report  investigates  the  problem  of  keeping  a  satellite  in  an 
orbit  around  the  earth  at  translunar  distances  from  the  earth.  The 
study  concentrates  on  distances  between  300*000  and  500,000  n  mi,  since 
in  this  range  the  earth’s  gravitational  attraction  is  the  dominant  fac¬ 
tor  in  producing  acceleration  of  a  satellite  vehicle  rela,xve  to  the 
earth.  The  approach  to  the  problem  involves  formulation  of  the  general 
nonlinear  equations  of  motion  for  an  object  under  the  gravitational  at¬ 
traction  of  the  earth,  the  sun,  and  the  moon.  The  solution  of  these 
equations  is  in  terms  of  the  variation  of  the  resulting  motion  from  a 
nominal  unperturbed  circular  orbit,  and  is  obtained  by  numerical  inte¬ 
gration  of  the  equations  of  motion.  The  results  indicate  no  instabil¬ 
ity  in  the  motion  of  the  orbital  plane,  which  is  characterized  by  a 
slow  regression  of  the  line  of  nodes  in  the  ecliptic,  while  the  inclin¬ 
ation  to  the  ecliptic  remains  essentially  constant.  The  investigation 
of  the  motion  in  the  orbital  plane  is  restricted  to  orbits  perpendicu¬ 
lar  to  the  ecliptic,  since  this  tends  to  reduce  the  likelihood  and  du¬ 
ration  of  near  approaches  of  the  satellite  to  the  moon,  thus  decreasing 
the  resulting  gravitational  disturbance.  Nevertheless,  the  in-plane 
motion  appears  to  be  characterized  by  an  instability  that  develops  as 
a  divergent  oscillation  in  the  radial  displacement  of  the  satellite 
from  the  reference  orbit.  The  onset  of  this  oscillation  is  found  to 
be  dependent  on  the  initial  geometry  of  sun,  moon,  and  satellite  rela¬ 
tive  to  the  earth  as  well  as  on  the  orbital  injection  conditions.  A 
method  is  developed  for  computing  these  injection  conditions  in  terms 
of  the  initial  geometry  so  that  the  development  of  instability  is  de¬ 
layed  for  two  to  three  years. 

Although  cases  are  found  in  which  the  radial  variations  from  a 
reference  orbital  radius  of  300,000  n  mi  remain  less  than  50,000  n  mi 
for  as  long  as  five  years,  it  appears  to  be  desirable  to  have  an  or¬ 
bital  control  capability.  This  is  done  by  a  modification  of  the  method 
of  determining  orbital  injection  conditions.  By  means  of  this  method 
a  satellite  can  he  kept  within  a  radial  displacement  of  26,000  n  mi  of 
the  reference  orbit  for  a  period  of  five  years  at  an  expenditure  of 
less  than  70  ft/sec  per  year  in  velocity  correction. 
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I .  INTRODUCTION 

If  the  total  acceleration  of  a  satellite  vehicle  relative  to  the 
earth  is  due  to  a  mutual  inverse  square  law  gravitational  attraction 
between  the  two  bodies,  the  motion  of  Hie  satellite  relative  to  the 
earth  is  said  to  be  Keplerian.  The  resulting  path  is  an  ellipse  with 
one  focus  at  the  earth's  center  and  a  fixed  orientation  relative  to 
inertial  space.  As  the  satellite  traverses  this  path,  its  orbital  an¬ 
gular  rate  varies  inversely  as  the  square  of  its  distance  from  the 
earth’s  center,  thus  maintaining  a  constant  orbital  angular  momentum. 

An  actual  earth  satellite  deviates  from  this  idealized  Keplerian 
motion,  since  its  acceleration  relative  to  the  earth  is  due  not  only 
to  the  earth's  inverse  square  law  gravitational  attraction  but  also  to 
perturbing  accelerations  caused  by  residual  atmospheric  drag,  solar 
radiation  pressure,  the  nonspherical  mass  distribution  of  the  earth, 
and  the  gravitational  attraction  of  the  sun  and  the  moon.  For  alti¬ 
tudes  up  to  atcut  100,000  n  mi,  all  of  these  perturbing  accelerations 
are  small  compared  to  the  inverse  square  term  in  the  earth’s  attrac¬ 
tion;  however,  their  relative  importance  varies  considerably  with  al¬ 
titude.  While  residual  atmospheric  drag  can  result  in  orbital  decay 
and  reentry  for  low-altitude  satellites,  its  effect  can  be  neglected 
at  altitudes  above  300  n  mi,  The  magnitude  of  the  acceleration  due  to 
solar  radiation  pressure  is  invariant  with  altitude,  depending  only  on 
the  area-to-mass  ratio  of  the  vehicle.  For  law  values  of  this  ratio 
the  radiation  pressure  effect  can  also  be  neglected. 

The  acceleration  due  to  the  nonspherical  mass  distribution  of  the 
earth  varies  inversely  as  the  fourth  or  higher  powers  of  the  distance 
from  the  earth's  center  and  is  negligible  at  distances  greater  than 
50,000  n  mi.  At  such  distances  from  the  earth  the  only  significant 
perturbing  accelerations  are  these  due  to  the  sur.  and  moon.  However, 
their  magnitudes,  even  at  a  distance  of  100,000  n  mi  from  the  earth, 
are  less  than  one  percent  of  the  earth’3  gravitational  attraction.  A 
perturbing  acceleration  of  this  magnitude  does  not  seem  likely  to  re¬ 
sult  in  orbital  instability,  and  an  examination  of  Kef,  1  indicates 
a  number  of  examples  of  satellites  at  altitudes  of  the  order  of 
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70,000  n  mi  that  have  been  in  orbit  for  8  to  12  years  with  no  apparent 
orbital  instability. 

As  the  orbital  radius  increases  and  becomes  comparable  with  that 
of  the  moon,  the  perturbing  acceleration  due  to  the  moon  may  become  as 
large  as  or  greater  than  the  earth's  gravitational  acceleration  at  times 
of  close  approach  of  the  satellite  to  the  moon.  Under  these  conditions 
the  vehicle  can  no  longer  be  regarded  as  an  earth  satellite.  However, 
it  can  be  shown  that  at  translunar  distances  of  the  order  of  300,000 
to  500,000  n  mi,  the  perturbing  accelerations  due  to  the  sun  and  the 
moon  are  again  small  compared  to  the  earth's  gravitational  attraction. 

This  report  presents  the  results  of  a  study  of  the  motion  of  a 
satellite  in  translunar  space  to  see  whether  a  stable  earth  orbit  can 
be  established  at  these  distances.  In  view  of  the  scarcity  of  experi¬ 
mental  data  on  translunar  earth-orbiting  vehicles,  the  study  is  of  ne¬ 
cessity  analytical  in  nature.  In  Section  II  the  general  equations  are 
developed  describing  the  motion  relative  to  the  earth  of  an  object  in 
translunar  space  under  the  Influence  of  the  gravitational  attractions 
of  the  earth,  the  sun,  and  the  moon.  In  Section  III  these  equations 
are  used  to  determine  the  variations  of  the  satellite  motion  from  a 
nominal  unperturbed  circular  orbit.  A  method  is  developed  for  the  de¬ 
termination  of  orbital  injection  conditions  that  minimize  these  varia¬ 
tions  and  produce  acceptable,  stable,  circular,  orbital  motion  for 
several  years.  In  the  event  that  the  variations  from  circularity  ex¬ 
ceed  permissible  tolerances,  a  method  is  also  developed  to  compute  the 
required  orbital  velocity  correction  and  its  time  of  application. 
Section  IV  presents  the  specific  conclusions  reached  in  this  study  in 
regard  to  the  feasibility  of  maintaining  tr»  tv*  lunar  orbits  either  pas¬ 
sively  or  actively,  while  some  of  the  details  of  the  analytical  devel¬ 
opment  of  the  orbital  injection  and  correction  technique  are  relegated 
to  Appendices  A  and  B.  The  computing  programs  used  are  Included  as 
Appendix  C,  together  with  a  definition  of  terms  and  a  description  of 
their  operation. 
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II.  METHOD  OF  ANALYSIS 

STATEMENT  OF  TOE  PROBLEM 

Tne  problem  investigated  in  this  report  can  be  stated  as  follows: 
Can  an  artificial  earth  satellite  he  established  in  a  stable  translunar 
earth  orbit  without  active  orbital  control,  and,  if  not,  what  is  the 
propulsion  requirement  to  maintain  such  an  orbit? 

iu  the  formulation  of  this  problem  certain  simplifications  and  as¬ 
sumptions  can  be  made  without  seriously  affecting  the  accuracy  of  the 
results.  The  acceleration  of  the  satellite  relative  to  the  earth  can 
be  expressed  as  the  vector  sura  of  the  following  elements:  (1)  the  ac¬ 
celeration  of  the  satellite  due  to  the  earth's  gravitational  attraction, 
which  can  be  represented  as  a  simple  inverse  square  law  variation  at 
the  assumed  translunar  distances  from  the  earth;  (2)  the  differential 
acceleration  of  the  satellite  and  the  earth  caused  by  the  moon's  gravi¬ 
tational  attraction;  and  (3)  the  differential  acceleration  of  the  sa¬ 
tellite  and  the  earth  caused  by  the  sun's  gravitational  attraction. 

The  motion  of  the  earth  relative  to  the  sun  is  assumed  to  be  an 

elliptical  orbit  in  the  plane  of  the  ecliptic  with  an  eccentricity,  e^., 

of  0.01674,  a  mean  orbital  angular  rate,  9Q,  of  0.0172  rad/solar  day, 

and  a  longitude  of  solar  perigee,  0p,  of  281*33' 53".  Similarly,  the 

motion  of  the  moon  relative  to  the  earth  is  also  assumed  to  be  an  el- 

* 

llptical  orbit  in  the  plane  of  the  ecliptic  with  an  eccentricity,  em, 
of  0.0549,  a  mean  angular  rate,  0 _ ,  of  0.22998  rad/solar  dav,  and  an 

©  iHO 

advance  of  perigee,  9  ,  of  0.001944  rad/solar  day. 

mp 

It  should  be  noted  that  in  this  representation  the  initial  posi¬ 
tions  of  the  earth  and  the  moon  in  their  orbits,  as  well  as  the  initial 
position  of  lunar  perigee,  are  independently  adjustable  rather  than  be¬ 
ing  Interrelated  through  an  ephemeris.  This  facilitates  the  study  of 
the  effect  of  the  initial  geometry  of  the  sun,  moon,  and  earth  on  the 
satellite  so  that  the  results  need  not  depend  on  the  particular  geome¬ 
try  that  may  exist  on  a  specific  launch  date. 

*lhis  Ignores  the  5*33'  inclination  of  the  moon's  orbit  to  the 
plane  of  the  ecliptic,  since  its  effect  on  the  satellite  orbit  is  small 
enough  that  the  additional  computation  Involved  did  not  seem  warranted. 
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COORDINATE  SYSTEMS 

In  Fig.  1  the  vector  relationships  of  the  sun,  earth,  moon,  and 
satellite  are  shown  schematically.  As  indicated  previously,  the  sun, 
moon,  and  earth  are  assumed  to  lie  in  the  plane  of  the  ecliptic;  how¬ 
ever,  there  is  no  such  restriction  on  the  satellite. 

The  coordinate  systems  used  in  the  formulation  of  the  equations 

of  motion  are  shown  in  Fig.  2  where  the  x  ,y  ,z  axes  define  an  earth- 

o  o  o 

centered  coordinate  system  with  no  rotation  relative  to  inertial  space. 

The  xoyo  plane  is  the  plane  of  the  ecliptic  and  the  x0  axis  is  in  the 

direction  of  the  first  point  of  Aries.  In  this  coordinate  system  the 

sun  and  moon  are  in  the  direction  of  the  lines  OS  and  OM  in  the  x  y 

o7o 

plane  making  angles  of  ©  and  6m  with  the  xq  axis.  The  x,y,z  axes  define 

an  earth-centered  coordinate  system  which  moves  in  such  a  way  that  Its 

x  axis  is  always  along  the  radius  vector  r  between  the  earth's  center 

and  the  satellite,  while  the  xy  plane  is  the  instantaneous  orbit  plane. 

The  orientation  of  this  system  relative  to  the  x  ,y  ,z  system  is  speei- 

o  o  o 

fied  by  three  angles,  T,  a,  and  6.  The  angle  ¥  is  the  orbital  regres¬ 
sion  angle  measured  from  the  x  axis  to  the  line  of  nodes  ON  at  the 

o 

Intersection  of  the  xy  and  xQyo  planes,  while  a  is  the  orbital  inclina¬ 
tion  angle  between  the  xy  and  xQyo  planes  and  0  is  the  angle  between 
the  x  axis  and  ON  measured  in  the  xy  plane. 

The  nine  direction  cosines  relating  the  unit  vectors  i,j,fc  of  the 
orbital  system  to  the  unit  vectors  Iq, jQ,ko  of  the  nonrotating  system 
are  given  in  Table  1  in  'terms  of  the  angles  V,  a,  and  0. 

In  addition,  it  is  convenient  to  define  two  unit  vectors  r,  and  I 

^  X  Ml 

in  the  direction  of  R  and  p,  respectively.  These  can  be  expressed  as 


?,  ■  I  cos  0  +  j  sin  0 
1  0  o 


I  ■  I  cos  9  +1  .*'in  0 
mo  m  Jo  m 


where  the  angles  0  and  fc  are  the  solar  and  lunar  longitudes. 


It  is  assumed  that  R  and  are  parallel. 


,s**?^!'***5,y*  ■*■ 


*<- 
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Table  1 

DIRECTION  COSINES 


The  acceleration  of  the  satellite  relative  to  the  earth  can  be 
expressed  in  vector  form  using  the  approach  of  Ref.  2  as  follows: 


tdicre  G  •  universal  gravitational  constant 
^  «  mss  of  the  earth 
Mh  »  mss  of  the  soon 
»  -  «g  ♦  1M/M, 

«£  »  semina jor  axis  of  the  earth’s  orbit 
aa  ■  tenlmajcr  axis  of  the  noon's  orbit 

The  principal  difference  between  Eq.  (3)  and  the  corresponding 
result  given  in  Ref.  2  as  Eq.  (B-37)  is  that  the  last  tern  in  Eq.  (3) 
is  an  exact  expression  for  the  differential  acceleration  of  the  satel¬ 
lite  and  tiie  earth  due  to  the  noon,  since  the  linearisation  of  this 


7  3v' 
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term  used  in  Ref.  2  is  not  valid  for  the  geometry  of  the  present 

problem.  However,  the  linearization  of  the  corresponding  solar  term 

is  still  valid  even  for  translunar  orbits. 

In  addition,  the  inclusion  of  the  eccentricity  of  the  orbits  of 

the  earth  and  the  moon  introduces  the  factor?#  ae/R  and  a _/p  which  can 

tt  n 

be  determined  as  a  function  of  time  by  the  following  equations: 


jj~  -  1  +  tE  cos  (0  -  9p) 


where 


”  “  1  +  e_  cos  (8  -  0  ) 

p  n  n  ip 


e. 


e  - 


b 


('  -  4) 


2W2  L1  +  ce  cog  (d 


-v]: 


0  * 


1  +  em  cos 


<»,- v] 


and  0  _  being  the  longitudes  of  solar  and  lunar  perigee. 
-  *P 

In  addition  the  ratio  p/r^  is  given  by  the  relation 


k-bf 


-1/2 


Hie  satellite  acceleration  can  also  be  expressed  as 

r(5  •S£)^li- 


d2r  j.2 

’  J?  '  *<VJ>  " 


+  \r?t  tr2<V£)]  +r<VI>(V5>(5 
+  7  ft  tr2(VJ)3  +  r(Sc*I)(;o*£)|C 


(4) 

(5) 


(6) 

(7) 


(8) 


(8) 


mwMW&W'yiWm *iwwws9!w^^  *  .-ro'^-'-  zr^rwrgW' 


./JWV^Vl-  >v*^v^.I'>-  vs«4  **»_»»  t^-J/  ^ 


—  >****,*  s’  •rC->-» 


niirfir  irri 


where  «o  is  the  angular  velocity  of  the  x,y,z  eye tea  relative  to  in¬ 
ertial  space ..  The  components  of  aQ  that  appear  in  Eq.  (9)  are  as 
follows: 


(a  •!)  ■  a  cos  G  +  y  sin  a  sin  e 
o 


(a  •?)  “  -a  sin  0  +  f  sin  a  cos  8 
o 


(a  »k>  ■  0  +  T  coe  « 
o 


However,  it  can  be  shown  that  if  the  xy  plane  is  the  instantaneous 
orbital  plane,  then  (wo»J)  aust  be  identically  zero.  Thus,  Eq.  (9) 


UMCO  • 

?  •  j^f  -  «vE>2)  +  7  h  p«vE>] 


j  +  r(5  -i)(w  -k)k  (13) 
o  o 


Combination  of  Eqs.  (10)  and  (11)  with  the  three  equations  obtained 
by  equating  components  of  tq*.  (3)  and  (13)  results  in  the  following 
equations  of  notion  for  the  satellite 


(r)3|+(l-^<vI)] 


h  <'V  -  ’^(r)5«! *  -r^irfi1  -  $)r  «V« 


T^ain.-sina  39jlr 


j^(r)3<V»<vs>  -  ¥(r)’(l  -  4)^  <vE>  j 
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Variatlonal  Equations 

A  considerable  improvement  in  computational  accuracy  can  be 
achieved  if  the  equations  of  motion  are  expressed  in  terms  of  the 
variations  of  r,  wz>  0,  ¥,  and  a  from  their  corresponding  values 
along  an  unperturbed  Keplerian  reference  orbit.  Since  there  is  no 
apparent  advantage  in  selecting  an  eccentric  reference  orbit,  it 
can  be  specified  as  follows: 


*  r 

0 

(19) 

• 

"  9o 

(20) 

• 

*  0 

0 

(21) 

-  0 

(22) 

«  0 

(23) 

where  the  mean  orbital  angular  rate, 


0Q,  is  given  by 


(24) 


In  addition,  it  is  also  convenient  to  use  0q  rather  than  the  time, 
t,  as  the  independent  variable  where  0q  is  given  by 


0 

o 


0  t 
o 


p 


(25) 
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and  represents  the  change  In  central  angle  along  the  reference  orhit 
during  the  time,  t. 

The  dependent  variables  r,  <oz,  6,  ¥,  and  a,  can  oe  expressed  in 
terms  of  their  variations  from  their  reference  orbit  values  by  the 
following  relations: 


r  ■  r  +  6r  *  r 
o  o 


H) 


a)  *9  +•$<0  *011  + 

z  o  *  o\ 


9*0  +60-0  /l  +  4~ - 

o  o\  d0 

\  o- 


¥  -  6¥  -  0 


a  ■  6a  *  0  "ts— 
o  d9 

o 


0  -  0  +  0(0)  +  60 
o 


¥  -  ¥(0)  +  5¥ 


a  *  a(0)  +  6a 


whe.--  y(0)4  ¥(0),  and  a(0)  are  the  initial  values  of  0,  ¥,  and  a 

when  ■'  (and  0  )  are  zero, 
o 

Similarly  the  angles  0,  0m>  and  0  can  be  expressed  in  terms 
of  the  independent  variable  0Q  by  the  relations 


®  “  .  6 o  +  «K°)  +  50 

0 

O 


(34) 


•(-*  -  ^ v  :^*ymnw*w»;s sss 
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0 

0  -  —  0  +  6  (0)  +  60m  (35) 

m  •  o  m  m 


0 

Q  =  JHE.  e+0  (0)  (3fi) 

mp  q  o  mp 

o 

where  0(0),  0^(0) ,  and  6mp(0)  are  the  initial  values  of  0,  0Jn,  and 
0  .  The  quantities  60  and  50m  are  the  angular  variations  of  the  sun 

and  the  moon  from  constant  angular  rate  motion  relative  to  the  earth. 

Substitution  of  Eqs.  (26)  through  (36'  into  Eqs.  (14)  through 
(18)  results  in  the  following  variational  form  of  the  equations  of 
motion. 


«*4*V 
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d6V  sin  9 


(40) 


30  ,  /a  3  _ 

jr(r)  ^i-ix^-k) 


y0 


3  v 

(r)( 

*-y 

rJ 

if  <VE> 

(41) 


where  the  dot  products  are  given  by  the  expressions 

"  cos  0  cos  (V  ~  0m)  _  c°s  «  sin  8  sin  (y  -  8  )  (42) 

jn 

*  “  sln  0  cos  (^  ~  em)  ~  cos  a  cos  8  sin  (y  -  8  )  (43) 


(1_ *k)  -  sin  a  sin  (y  -  8  ) 
ra  m 

(^l*i)  "  cos  8  cos  (V  -  0)  -  cos  a  sin  8  3in  (y  -  0) 
(rj*J)  “  “  sin  8  cos  (y  -  0)  -  cos  a  cos  8  sin  (y  -  0) 


(44) 

(45) 

(46) 


(r^-k)  -  sin  a  sin  (y  -  0) 


(47) 


Finally,  the  values  of  60  and  68m  can  be  obtained  from  Eqs.  (6) 
and  (7)  after  substitution  of  Eqs.  (34)  through  (36)  to  give 


w!  jr  r-^  $-<&*»  r»$ 
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d60 

d6 

o 


dse 

_ m 

d0 

o 


(48) 


(49) 


Equations  (37)  through  (41)  represent  the  complete  equations  of 
motion  applica’  le  to  satellites  in  translunar  orbits.  They  consti¬ 
tute  a  set  of  five  coupled  nonlinear  differential  equations  in  the 

variables  6r,  6u  ,  60,  df,  and  6a  as  functions  of  the  independent 
z 

variable  0  .  It  should  be  emphasized  that  this  formulation  puts  no 
o 

limitation  on  the  magnitude  of  any  of  the  five  dependent  variables. 

Determination  of  the  orbital  motion,  in  view  of  the  nonlinear¬ 
ity  of  the  equations,  requires  simultaneous  numerical  integration 
not  only  of  the  five  equations  of  motion  but  also  of  Eqs.  (48)  and 
(49)  to  determine  the  values  of  60  and  60^  as  a  function  of  the  in¬ 
dependent  variable  0Q.  While  this  can  be  done,  it  is  found  that  cer¬ 
tain  approximations  can  be  made  that  markedly  reduce  the  solution 
time  without  any  significant  change  in  the  nature  of  the  results. 


Uncoupled  Equations 

Equations  (40)  and  (41)  describe  the  part  of  the  motion  that  re¬ 
sults  in  a  change  in  the  orientation  of  the  orbital  plane,  while  Eqs. 
(37)  through  (39)  determine  the  motion  in  the  orbital  plane  that  re¬ 
sults  in  changes  of  size  and  shape  of  the  orbit.  While  these  two 
types  of  motion  are  not  completely  uncoupled,  the  coupling  is  weak 
enough  so  that  the  motion  of  the  orbital  plane  can  be  treated  sepa¬ 
rately  by  ignoring  the  effects  of  6r,  60,  6u>  ,  and  6a  on  the  right 
side  of  Eqs.  (40)  and  (41).  With  the  additional  assumption  that  em 

and  are  zero,  these  equations  become 

E 
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(51) 

These  equations  can  be  integrated  numerically  along  with  Eqs.  (48) 
and  (49)  to  determine  the  variations  in  the  orientation  of  the  orbital 
plane  as  given  by  6'P  and  5a.  It.  is  found,  as  will  be  shown  in  the  sec¬ 
tion  on  results,  that  the  solution  for  541  consists  of  a  secular  term 
with  small  oscillatory  variations  superposed,  while  6a  can  be  repre¬ 
sented  by  a  small  negative  bias  term  with  small  oscillatory  terms 
superposed. 

In  a  similar  manner,  the  variation  of  the  in-plane  motion  can 

also  be  determined  independently  by  replacing  a  with  its  reference 

value,  a(0),  and  V  with  its  secular  variation  in  Eqs.  (37)  through  (39). 

These  equations  can  also  be  integrated  numerically  along  with  Eqs.  (48) 

and  (49)  to  give  the  in-plane  variations  of  the  motion  as  specified 

by  6r,  60,  and  6<jj  . 

z 

The  in-plane  motion  can  be  further  simplified  for  two  special 
cases  that  are  of  interest,  namely  when  the  orbital  inclination  to 
the  ecliptic  is  either  0  or  90  deg. 

If  a(0)  is  equal  to  0  deg,  the  orbit  lies  in  the  plane  of  the  eclip' 
tic  and  the  line  of  nodes,  ON,  in  Fig.  2  is  indeterminate.  As  a  re¬ 
sult,  the  angles  8  and  as  well  as  their  variations  60  and  AY,  are 
also  indeterminate.  However,  the  sum  of  0  and  ¥  is  determinate  and 
represents  the  angular  displacement  of  the  satellite  about  the  earth's 

center  measured  from  the  x  axis.  If  his  angle  is  defined  as  y,  then 
•  " 

by  Eq.  (18)  y  is  equal  to  u  and  Eq.  (39)  reduces  to  an  identity, 

z 

while  Eqs.  (37)  and  (38)  become 
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Thus  the  ln-plane  motion  is  completely  described  by  solving 
these  equations  for  5r  and  since  0  and  ¥  no  longer  occur  sepa¬ 
rately  in  the  equations. 

If  a  is  90  deg,  then  by  Eq.  (18)  0  is  equal  to  u  ,  and  Eq.  (39) 
o  « 

again  reduces  to  an  identity,  while  Eqs.  (37)  and  (38)  take  the  form 


■p 
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9*y 

3©:  /a„\3 


L  +  djeN  d  (6r\  _  ^(!e\  sln  20  cog2  Q 

(l  +  k j  V  dV  deo  VV  20^  Vr  ' 


!«■)’(■  -4)| 


sin  0  cos  0 


Thus  the  in-plane  motion  can  again  be  described  by  solving  Eqs. 
(54)  and  (55)  for  the  variations  6r  and  60. 


i  %?"•>» -4 
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III.  RESULTS  ANI>  DISCUSSION 
SELECTION  OF  A  REFERENCE  ORBIT 

If  the  orbit  of  .  satellite  .bout  the  e.rth  i.  to  be  stable,  “ 
1,  necessary  that  Us  acceler.tion  due  to  the  earth's  5r.eit.tt0u. 
attraction  be  large  conb.red  ulth  the  differential  acceleration.  o 
the  satellite  end  the  earth  caused  by  either  the  sun  or  the  aoon. 
Thus  an  euaminatien  of  the  masnitudes  of  these  accelerations,  a.  a 
function  of  distance  fro.  the  earth,  should  give 
selection  of  the  orbital  radius  of  the  reference  orbit 
the  orbital  variations  are  to  be  neaaurad.  The  eaxims.  differential 
acceleration  of  the  satellite  end  the  earth  due  to  the  -ecu  occur, 
when  the  satellite  end  the  noon  have  a  co-on  faction  froe  he 

earth's  center  (I  -  i„>.  *n  —^tlon  »'  <»  **"  * 

oaxinum  acceleration  1.  along  t  and  ha.  a  nagnitude.  A,,,  give 


e‘ 


wr 


-  r 


*0 

w 


(‘•I- 


(56) 


.here  the  eccentricity  of  the  -con's  orbit  !•  neglected. 

The  acceleration  of  the  satellite  by  the  earth,  V  i»  **«•»  • 
the  expression 


**"T 


(57) 


Thus  the  relative  Inherence  of  the  lunar  effect  can  be  detemlned 
as  the  ratio  of  A^  to  Ag  as  follows: 


(58) 


Since  Is  given  by  the  relation 

InC 
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Eq.  (58)  becomes 


f2 

V;  ^  r  (r  ~  S> 

*E  '  «E  r3 
U  m 


(59) 


(60) 


For  1  equal  to  i  ,  r  is  equal  to  the  absolute  value  cf  (r  -  a  )  and 

ah  a 

Eq.  (60)  can  be  written  as 


2 

r 


(r  - 


(6.1) 


where  the  plus  sign  applies  if  r  is  greater  than  aR  and  the  minus  sign 
for  r  less  than  a  . 

A 

Similarly,  the  maximum  differential  acceleration  of  the  satellite 
and  the  earth  due  to  the  sun  occurs  when  the  directions  of  the  satel¬ 
lite  and  the  sun,  relative  to  the  center  of  the  earth,  are  either  the 
same  or  exactly  opposite  (t  *  >r  j) .  In  either  case  Eq.  (3)  shows  that 
this  maximum  acceleration  is  along  r  and  has  a  magnitude,  Agf,  given  by 

Sr  "  2*o  (62) 


where  the  eccentricity  of  the  earth's  orbit  is  also  neglected.  Since 
8o  is  given  by  the  relation 


the  ratio  of  Ag  to  Ag  can  be  expressed  in  the  form 


Sr  f!Wr\3 

*E  '  *E  W 


(63) 


(64) 
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only  during  close  approaches  of  the  satellite  to  the  moon*  During  the 
rest  of  the  orbit  the  average  lunar  effect  is  of  the  order  of  6  percent 
or  about  twice  the  average  solar  effect.  Thus  if  the  frequency  of  oc¬ 
currence  and  the  duration  of  these  close  approaches  of  the  satellite 
to  the  noon  can  oe  appreciably  reduced  b?  a  suitable  choice  of  the  or¬ 
bital  orientation,  then  the  average  disturbance  due  to  the  noon  can  be 
reduced  and  the  possibility  of  a  stable  orbit  enhanced. 

Thus  the  refetunce  orbit  used  for  most  of' this  study  is  a  circular 
orbit  with  a  radius  of  300,000  n  al  and  a  period  of  47.8  days.  The 
variation  of  this  orbital  period  with  orbital  radius  In  the  translunar 
region  is  shown  in  Tig.  4. 


COMPUTATIONAL  PKOCEPCTE 

In  the  previous  section, 
it  was  shown  that  a  range  of 
t ran s lunar  distances  exists  in 
which  the  lunar  and  solar  ef¬ 
fects  are  win  1*1  zed  relative 
to  the  earth's  attraction. 
However,  it  is  still  necessary 
to  de  re  twine  whether,  wider 
the  influence  of  these  resid¬ 
ual  disturbances,  a  suitable 
stable  oxblt  can  exist.  This 
involves  the  solution  of  the 
equations  of  notion  derived 
in  Section  II.  In  view  of  the 
nonlinearity  of  these  equations, 
it  is  obvious  that  the  solution 
must  he  obtained  by  numerical 

Integration.  This  is  accomplished  by  programing  the  equations  for 

* 

solution  on  the  JOSS  computer  using  a  fourth  order  RungeHCutt* 

*J0SS  Is  the  trademark  of  the  remote  console,  tine  sharing  com¬ 
puter  -ystem  developed  at  The  Hand  Corporation. 


Fig.  4 — Orfciu!  pcjtkxl  xi  ?  feactvn  of 
orbital  radios 
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integration  procedure.  Initially,  the  completely  general  equations 
as  represented  by  Eqs  (37)  through  (41)  were  programmed.  However, 
this  involves  the  simultaneous  integration  of  these  five  equations  as 
well  as  Eqs.  (43)  and  (49)  for  the  motion  of  the  3un  and  the  moon,  and 
the  solution  time  is  found  to  be  excessive.  If  the  uncoupled  form  of 
the  equations  of  motion  is  used,  then  the  determination  of  the  motion 
of  the  orbital  plane  involves  the  integration  of  Eqs.  (50)  and  (51)  to 
obtain  64*  and  6a,  as  well  as  Eqs.  (48)  and  (49)  to  obtain  60  and  69^. 
The  in-plane  motion  is  determined  independently  by  integrating  the  un¬ 
coupled  forms  of  Eqs.  (37)  through  (39)  to  obtain  6r,  6w  ,  and  60. 

z 

Comparison  of  the  lesults  obtained  from  the  uncoupled  equations  with 
those  from  the  more  general  form  shows  no  significant  difference  in 
the  nature  of  the  orbital  motion.  Thus,  the  uncoupled  form  of  the 
equations  is  used  in  all  cases. 

MOTION  OF  THE  ORBITAL  PLANE 


Motion  Relative  to  the  Ecliptic  Plane 

An  examination  of  Eqs.  (50)  and  (51)  shows  that  the  right-hand 
side  of  Eq.  (50)  averaged  over  0q  has  a  nonzero  value,  while  a  similar 
average  of  the  right  side  of  Eq.  (51)  is  zero.  Thus,  when  the  two 
equations  are  Integrated,  641  exhibits  a  secular  variation  with  small 
oscillatory  variations  superposed,  while  6a  undergoes  small  oscilla¬ 
tions  about  a  small  negative  bias.  Thus,  to  a  good  approximation,  the 
motion  of  the  orbital  plane  can  be  described  as  a  slow  regression  of 
its  line  of  nodes  in  the  plane  of  the  ecliptic,  while  its  inclination 
to  the  ecliptic  remains  constant.  Figure  5  is  a  plot  of  the  regres¬ 
sion  period  as  a  function  of  orbital  inclination  to  the  ecliptic  for 
a  300,000-n  mi  orbital  radius.  It  is  seen  that  the  period  increases 
monotonically  with  inclination  heving  a  value  of  about  five  years  at 
0  deg  and  becoming  infinite  at  90-deg  Inclination  to  the  ecliptic. 

Thus  for  a  90-deg  inclination  the  line  of  nodes  simply  executes  small 
angular  oscillations  about  a  fixed  direction  in  the  plane  of  the  eclip¬ 
tic.  To  give  some  idea  of  the  variation  of  regression  period  with 
altitude,  several  points  are  shown  in  Fig.  5  for  an  orbital  radius 
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of  400,000  n  mi.  It  is  seen 
that  for  a  given  inclination 
the  regression  period  decreases 
with  altitude,  particularly  at 
higher  values  of  orbital  in¬ 
clination. 

In  Fig.  6  a  plot  of  the 
bias  of  the  orbital  inclina¬ 
tion  angle  is  shown  as  a  func¬ 
tion  of  the  nominal  Inclination, 
a(G).  The  bias  also  varies  mono- 
tonically  from  a  value  of  '  deg 
when  a(0)  is  zero,  on  t  a  value 
of  about  -2.5  deg  when  a(0)  is 
equal  to  90  deg.  As  in  Fig,  5, 
several  points  are  shown  for 
an  orbital  radius  of  400,000 
n  mi  and  it  is  seen  that  the 
magnitude  of  the  bias  increases 
with  altitude.  Since  the  am¬ 
plitudes  of  the  oscillatory  terms  in  both  SV  and  6a,  as  well  as  the 
steady-state  regression  rate  appear  to  be  bounded,  the  resulting  mo¬ 
tion  of  the  orbital  plane  does 
not  constitute  an  orbital 
instability. 


Motion  Relative  to  the 
Equatorial  Plane 

In  the  case  of  earth  sat¬ 
ellites,  the  earth's  equatorial 
0  iojo3o«omwa)»o9o  plane,  rather  than  the  ecliptic, 

Ofbitol  inciinotion  to  th«  elliptic  { d*g > 

is  ordinarily  used  as  a  refer¬ 
ence  relative  to  which  the  mo¬ 
tion  of  the  orbital  plane  is 
determined.  If  this  is  done  in 


Fig.  5 — Regression  period  as  a  function  of 
orbital  inclination 


Fig.  fi — Orbital  inclination  bias  as  a 
function  of  orbital  inclination 
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the  case  of  the  motion  described  in  the  previous  section,  it  is  found 
that  neither  the  orbital  inclination  nor  the  regression  rate  of  the 
node  remain  constant  in  this  equatorial  reference  system.  As  an  example 
an  orbit  that  is  initially  in  the  equatorial  plane  is  inclined  at  an 
angle  of  23*27'  to  the  ecliptic.  From  Fig.  5  the  regression  period  of 
the  node  in  the  ecliptic  is  6.7  years;  thus  after  only  3.35  years,  th? 
orbital  inclination  relative  to  the  equatorial  plane  would  increase 
from  0  deg  to  46 °54/ .  The  orbital  inclination  to  equatorial  plane 
can  remain  constant  only  if  its  orientation  relative  to  Inertial  space 
remains  constant.  This  condition  is  satisfied  if  the  orbital  inclina¬ 
tion  to  the  ecliptic  is  either  0  or  90  deg.  If  the  inclination  to  the 
ecliptic  is  0  deg,  the  corresponding  Inclination  to  the  equator  Is 
uniquely  determined  as  23*27' .  On  the  other  hand,  if  the  inclination 
to  the  ecliptic  is  90  deg,  then  the  corresponding  inclination  to  the 
equatorial  plane,  a_,  is  given  by  the  relation 

\7 

cos  a„  -  cos  ¥(0)  sin  X  (65) 

It 

where  Y(0)  ■  steady-state  position  of  the  line  of  nodes  In  the  ecliptic 
X  ■  angle  between  the  equatorial  plane  and  the  ecliptic. 

Thus,  with  a  90-deg  inclination  to  the  ecliptic,  the  constant  equato¬ 
rial  inclination  may  have  a  value  anywhere  from  66*33'  up  to  a  90-deg 
polar  orbit. 

Although  it  may  not  be  necessary  to  hold  a  fixed  orbital  inclina¬ 
tion  to  the  equator,  it  does  appear  that  an  inclination  cf  90  deg  to 
the  ecliptic  has  certain  advantages  in  minimizing  the  variations  of 
the  in-plane  motion  resulting  from  close  approaches  of  the  satellite 
to  the  moon.  This  is  discussed  further  in  the  next  section. 

Preferred  Orbital  Inclination 

If  a  satellite  is  in  a  nominal  300,000-n  mi  orbit,  Fig.  3  shows 
that  the  maximum  perturbing  acceleration  due  to  the  moon  is  15.  per¬ 
cent  of  the  acceleration  due  to  the  earth.  This  maximum  acceleration 
occurs  whenever  the  satellite  and  the  moon  are  in  the  same  direction 
relative  to  the  earth.  At  this  time  the  separation  of  the  satellite 
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and  the  moon  is  equal  to  the  difference  of  their  orbital  radii,  92,572 
n  mi,  if  the  eccentricity  of  the  moon's  orbit  is  neglected.  If  the 
satellite's  orbit  and  the  moon’s  orbit  lie  in  the  same  plane,  then  one 
of  these  close  approaches  occurs  each  time  the  moon  gains  one  complete 
revolution  relative  to  the  satellite.  This  occurs  with  a  periodicity., 
of  63.47  days  corresponding  to  the  difference  of  the  two  orbital  angu¬ 
lar  rates.  The  variation  of  the  radial  component  of  the  moon's  per¬ 
turbing  acceleration  during  one  of  these  close  approaches  is  shown  in 
Fig.  7,  which  is  a  plot  of  the  lunar  term  of  Eq.  (52)  as  a  function  of 
time,  where  zero  time  occurs  at  closest  approach. 

If  the  satellite's  orbital  plane  has  an  inclination  of  90  deg  to 
that  of  the  moon,  a  close  approach  to  the  moon  can  only  occur  at  the 
two  points  of  intersection  of  the  satellite  orbit  with  the  plane  of 
the  moon's  orbit.  If  the  orbital  periods  of  the  moon  and  the  satel¬ 
lite  are  incommensurate,  the  probability  of  a  close  approach  is  much 


T»*»*  (rfvyi) 


Fig.  7  Fig.  8 

Radial  component  of  lunar  perturbing  acceleration  as  a  function  of  time 


less  than  in  the  previous  case  where  it  could  occur  at  any  point  along 
the  satellite  orbit. 


However,  if  both  the  satellite  and  the  moon  do  arrive  at  the  line 
of  intersection  of  the  two  orbital  planes  simultaneously,  then  the  vari¬ 
ation  of  the  radial  component  of  the  perturbing  acceleration  with  time 
is  shown  in  Fig,  8,  which  is  a  plot  of  the  lunar  term  in  Eq.  (54).  It 
is  seen  that  although  the  peak  acceleration  ii  the  same  as  that  in  Fig. 

7,  the  width  of  the  peak  is  considerably  less.  Thus,  the  total  Impulse 
as  measured  by  the  area  under  the  curve  is  less.  In  addition.  Fig.  9 
shows  the  effect  on  this  acceleration  peak,  if  the  arrival  of  the  moon 
at  the  intersection  differs  from  that  of  the  satellite.  If  the  moon  ar¬ 
rives  earlier  than  the  satellite,  the  peak  occurs  slightly  before  the 
satellite  reaches  the  intersection  and  the  height  of  (he  peak  is  reduced. 
The  height  of  the  peak  is  similarly  reduced  for  a  late  arrival  of  the 
moon,  but  it  occurs  after  the  satellite  has  passed  the  intersection. 

On  the  basis  of  a  comparison 


Kig.  9 — Variation  of  lunar  acceleration 
with  difference-  in  arrival  time  (90* 
orbital  inclinati  uj) 


of  Figs.  7,  8,  and  9,  it  is 
evident  that  the  perturbing 
acceleration  due  to  the  moon 
is  less  severe  for  the  90-deg 
inclination.  Thus,  in  the  in¬ 
vestigation  of  the  in-plane 
motion  described  in  the  next 
section,  the  principal  et  *ha- 
sis  is  on  orbits  with  this  in¬ 
clination  to  the  ecliptic. 


MOTION  IN  THE  ORBITAL  PLANE 


Preliminary  Results 


The  undisturbed  in-plane 
motion  of  the  satellite  along 
the  300,000-n  mi  radius  refer¬ 
ence  orbit  has  an  orbital  ve¬ 


locity  of  2774.4  ft/sec,  which 
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results  In  an  orbital  angular  rate  of  0.13142  rad/day  and  an  orbital 
period  of  47.81  days.  The  orbital  inclination  of  90  deg  to  the  eclip¬ 
tic  precludes  any  significant  change  in  orientation  of  the  plane  as 
shown  previously.  The  actual  motion  is  determined  relative  to  the 
reference  satellite  by  numerical  Integration  of  Eqs.  (54)  and  (55) 
for  6r  and  59,  which  represent  the  radial  and  tangential  displacements 
of  the  actual  satellite  from  its  reference.  At  the  same  time,  inte¬ 
gration  of  Eqs.  (48)  and  (49)  determines  the  motion  of  the  sun  and 
the  moon  relative  to  the  earth.  The  input  data  required  include  the 

initial  values  of  0,  6,  0  ,  5r,  dSr/d0  ,  and  d60/d0  ,  the  initial  val- 

moo 

ues  of  50,  50,  and  50  being  equal  to  zero.  The  equations  were  inte¬ 
rn 

grated  with  a  0.5-day  Integration  step  size,  while  the  result  is 

printed  out  at  one-day  intervals.  The  basic  output  of  the  program 

gives,  as  a  function  of  time,  the  quantities  6r/rQ,  60,  and  their 

rates  of  change  with  respect  to  0Q  as  well  as  the  variations  60  and 

60  .  These  quantities,  together  with  characteristics  of  the  reference 
m 

orbit,  are  then  sufficient  to  completely  describe  the  motion  of  the 
actual  satellite  relative  to  the  earth.  The  description  of  the  pro¬ 
gram  is  given  in  Appendix  C, 

The  ln-plane  motion  of  a  satellite  can  be  regarded  as  stable  if 
the  magnitude  of  6r  never  exceeds  some  preassigned  value,  thus  main¬ 
taining  a  reasonable  approximation  to  a  circular  orbit.  However,  the 
cost  of  demonstrating  any  such  absolute  stability  cf  an  orbit  by  ma¬ 
chine  solution  becomes  exorbitant.  For  the  purpose  of  this  report 
the  in-plane  motion  is  consequently  regarded  as  stable  if  the  magni¬ 
tude  of  6r  remains  less  than  50,000  n  mi  for  a  period  of  five  years. 

As  a  first  step  in  the  study  of  this  stability,  it  is  of  interest 
to  determine  triiether  the  resulting  motion  is  sensitive  to  the  initial 
geometry  of  the  satellite,  the  sun,  and  the  moon  relative  to  the  earth 
as  specified  by  the  initial  values  of  the  angles  0,  0,  and  0  (,  This 
is  investigated  by  varying  the  initial  value  of  0  in  steps  of  10  deg 
for  the  same  initial  positions  of  the  sun  and  the  moon.  The  initial 

values  of  6r,  d6r/d0  ,  and  d69/d0  are  oet  equal  to  zero  with  the  re- 

o  o 

suit  that  in  all  cases  the  orbital  injection  conditions  are  those  of 
the  reference  orbit. 
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An  examination  of  the  resulting  orbital  motion  shows  a  signifi¬ 
cant  variation  depending  on  the  position  of  orbital  injection  as  speci¬ 
fied  by  the  initial  value  of  9.  In  Fig,  10  the  variation  of  the  radial 
deviation,  Sr,  from  the  reference  orbit  is  shown  as  a  function  of  time 
for  two  cases  in  which  the  only  difference  is  an  initial  value  of  190 
deg  for  G  in  the  upper  curve  and  a  value  of  20  deg  in  the  lower  curve. 
These  two  curves  are  roughly  the  extremes  of  the  types  of  notion  that 
result.  In  the  upper  curve  representing  the  worst  case,  an  oscillation 
in  Sr  develops,  reaching  in  amplitude  in  excess  of  50,000  n  mi  within 
250  days,  while  the  frequency  of  the  oscillation  is  approximately  that 
of  the  reference  orbit.  In  the  lower  curve  of  Fig.  10,  which  represents 


(a) 
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Time  (days) 

r0  =  300,000  n  mi 
9m  (0)  =  225° 
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Fig.  10 — Effect  of  initial  geometry  on  orbital  stability 
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the  best  case,  the  amplitude  of  the  oscillation  in  Sr  does  not  in¬ 
crease  as  rapidly  and  its  periodicity  is  more  complex,  involving  os¬ 
cillations  at  frequencies  higher  than  that  of  the  reference  orbit. 

In  addition  to  the  radial  variations  shown  in  Fig.  10  there  are 
also  tangential  variations,  which  in  both  cases  are  characterized  by 
a  slow  dirft  of  the  satellite  ahead  of  the  corresponding  position  on 
the  reference  orbit.  This  constitutes  a  steadv-state  increment  in 
the  orbital  ?.ngular  rate,  which  can  be  explained  qualitatively  by  the 
fact  that  the  average  radial  acceleration  of  the  satellite  relative 
to  the  earth  is  increased  by  the  presence  of  the  sun  t;»d  the  moon. 
However,  this  change  in  orbital  angular  rate  is  not  regarded  as  an 
orbital  Instability,  since  it  does  not  alter  the  shape  of  the  orbit. 
On  the  other  hand,  the  radial  variations  shewn  in  the  upper  curve  of 
Fig.  10  have  exceeded  the  50,000-n  mi  stability  limit  in  a  little 
over  160  days. 

If  a  value  of  -50,000  n  mi  for  6t  occurs  at  the  time  of  closest 
approach  of  the  satellite  and  the  noon,  the  satellite  is  about  42,000 
n  mi  from  the  moon  and  250,000  n  mi  from  the  earth.  An  examination 
of  Fig.  3  shows  that  under  these  conditions  the  lunar  gravitational 
disturbance  is  about  45  percent  of  that  due  to  the  earth  as  compared 
with  15  percent  at  300,000  n  mi.  Such  an  Increase  in  the  disturbing 
acceleration  can  result  in  even  greater  divergence  in  the  oscilla¬ 
tion  of  Sr. 

From  this  discussion  it  is  seen  that  the  resulting  orbital  mo¬ 
tion  is  a  function  of  the  initial  geometry  of  the  satellite,  the  sun, 
and  the  moon  relative  to  the  earth.  While  it  is  evident  that  insta¬ 
bility  can  occur,  it  is  not  necessarily  Inevitable.  Hopefully,  by  a 
suitable  choice  of  initial  geometry  and  orbital  Injection  conditions, 
a  long  duration  stable  orbit  may  be  possible. 

Long  Duration  Otblt 

If  the  solution  for  case  (b)  shown  in  Fig.  10  is  extended  beyond 
the  250  days,  the  resulting  radicl  variations  are  shown  in  the  posi¬ 
tive  time  region  in  Fig.  11.  It  le  again  seen  that  an  oscillatory 
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dlvergence  develops  at  orbital  frequency,  and  Its  amplitude  exceeds 
50, 000  n  *1  after  471  days.  Since  the  computer  has  no  compunction 
about  running  tine  backward,  it  is  possible,  by  using  a  negative  com¬ 
puting  step  site,  tc  determine  the  variation  of  6r  that  would  haw 
occurred  prior  to  zero  tine  in  order  to  achieve  the  specified  initial 
conditions  at  zero  time.  The  resulting  solution  is  shown  along  the 
negative  tine  axis  in  Fig.  11.  The  magnitude  of  6r  remains  less  than 
50,000  n  mi  for  1362  days  in  the  negative  time  direction.  The  values 
of  6r,  dSr/d®  ,  66,  dS6/dO  ,  60,  and  66  corresponding  to  any  point 

O  C  m 

on  this  curve  constitute  a  set  of  Initial  conditions  for  orbital  in¬ 
jection  on  that  day  which  will  generate  the  curve  in  Fig.  11  to  the 
right  of  the  Injection  point.  Thus,  if  a  satellite  had  been  Injected 
into  orbit  at  time  -1362  days  with  the  initial  conditions  correspond¬ 
ing  to  that  point  in  Fig.  11,  the  resulting  6r  variation  would  have  a 
magnitude  less  than  50,000  n  ml  until  +471  days,  a  total  Interval  of 
1833  days  or  5.02  years,  thus  satisfying  the  previous  definition  of 
orbital  stability.  This  procedure  of  determining  initial  conditions 
requires  an  accurate  epheaerls  of  the  motion  of  the  moon  and  the  sun 
relative  to  the  earth  for  at  least  five  years  in  advance  of  the  pro¬ 
posed  launch  date.  In  addition,  a  considerable  amount  of  computing 
time  might  be  required  to  find  a  time  in  the  future  from  which  a  stable 
orbit  could  be  established  over  a  period  including  the  five-year  in¬ 
terval  starting  from  the  desired  launch  time.  Finally,  even  after  the 
initial  conditions  are  determined  the  precision  required  in  realising 
these  conditions  so  that  the  orbit  would  follow  the  computer  version 
for  five  years  might  be  difficult  to  realize. 

Determination  of  Initial  Conditions 

It  is  found  that  by  making  small  changes  in  the  orbital  Injection 
conditions  from  those  corresponding  to  the  reference  orbit,  consider¬ 
able  Improvement  in  the  resulting  orbital  motion  can  be  achieved. 
However,  the  use  of  an  empirical  cut  and  try  method  to  improve  the 
injection  conditions  is  rather  uncertain  and  very  tine  consuming. 

While  the  specification  of  the  optimum  initial  conditions  for 
the  complete  nonlinear  equations  of  motion  is  very  difficult,  it  Is 
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possible  to  make  use  of  the  linearized  fora  of  these  equations  for  such 
a  determination.  If  the  higher  order  terns  in  Eqs0  (54)  and  (S5)  are 
neglected,  the  result  is  the  following  set  of  linearized  equations: 


where  the  eccentricity  of  the  orbits  of  the  earth  and  the  noon  are  ig¬ 
nored  .  In  addition,  it  is  assumed  €6  remain*  c^all  so  that  it  can  be 
neglected  on  the  right  side  of  Eqs.  (6$)  and  (67). 

Appendix  A  shown  that  the  terms  on  the  right  side  of  Eqs.  (66) 
and  (67)  can  be  expressed  by  a  series  expansion  as  follows: 
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where  Che  convergence  of  the  aeries  in  Eqs..  (68)  end  (61')  Is  excellent 

for  s  value  of  20  for  n  . 

o 

Substitution  for  6,  9,  end  0^  in  terns  of  0Q  by  means  of  Eqs.  (31) » 
(36),  and  (35)  transforms  Eqs.  (66)  and  (67)  to 
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Thu.s,  the  driving  functions  of  the  equations  of  motion  are  expressed 
as  summations  of  sinusoidal  terms  whose  arguments  are  functions  of  0  . 

By  superposition,  the  solution  of  Eqs.  (72)  and  (73)  can  be  ex¬ 
pressed  as  the  summation  of  the  solutions  for  the  individual  driving 
terms  in  the  form 
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(In  Eqs.  (76),  (77),  (78),  and  (81)  the  summations  do  not  include 
the  term  for  i  =  j  =  0  if  either  or  occurs  as  a  factor  in  the 

denominator.) 
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An  examination  of  the  linearized  solutions  for  fir/r  and  dfi0/d0 

o  o 

given  by  Eqs.  (76)  and  (77)  shows  that  the  choice  of  the  initial  con- 
*  /  / 

ditions  firQ,  6rQ,  and  50q  appear  in  the  constant  terms  and  the  coef¬ 
ficients  of  sin  0  and  cos  0  but  not  in  the  coefficients  of  the 
o  o 

driven  solution  a^,  b^,  c^,  and  d^. 

If  these  initial  conditions  are  set  equal  to  zero,  the  orbital 

injection  conditions  are  those  for  the  reference  orbit.  Under  these 

conditions,  fir^,  the  constant  bias  term  in  fir,  is  equal  to  r  A  , 
j  o  o 

while  B  is  the  steady-state  rate  term  in  d56/d0  .  The  amplitude  of 

O'  o 

the  orbital  frequency  oscillatory  term  in  fir/rQ  is  given  by 


5r 


osc 


V*l  +  A2 
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and  50 ' 
o 


are  the  initial  values  of  dfir/d0 
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and  d(50)/d0  . 
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while  the  corresponding  amplitude  in  d60/d0Q  is  twice  this  value. 

For  0(0)  equal  to  zero  the  value  of  this  oscillatory  amplitude  in 

6r/ro  due  to  the  lunar  terms  is  plotted  in  Fig.  12a  as  a  function  of 

the  initial  lunar  longitude  0  (0);  Fig.  12b  is  the  amplitude  due  to 

m 

the  solar  terms  plotted  as  a  function  of  the  initial  solar  longitude 
0(0).  The  curves  are  plotted  only  up  to  180  deg,  since  both  are  even 
functions  of  the  respective  latitudes.  The  reason  for  the  observed 
sensitivity  of  the  orbital  notion  to  the  initial  geometry  of  the  sun, 
moon,  and  satellite  is  evident  from  Figs.  12a  and  12b. 

However,  considerable  improvement  in  the  orbital  motion  can  be 
achieved  by  an  appropriate  choice  of  initial  conditions.  Since  the 
development  of  orbital  instability  seems  to  be  characterized  by  the 
buildup  of  a  divergent  oscillation  in  dr  at  orbital  freque.  'ey,  it 
seems  desirable  to  reduce  the  initial  amplitude  of  this  oscillation 
to  zero.  This  can  be  accomplished  in  both  Eqs.  (76)  and  (77)  by  sat¬ 
isfying  the  relations 


— -  +  260'  - 
ro 


dr' 


In  addition,  the  assumption  that  60  remains  small  cannot  be  satisfied 
if  there  is  a  steady-state  term  in  Eq.  (77).  This  can  also  be  made 
zero  by  satisfying  the  relation 


dr  B 

2  — —  +  d0'  -  ~ 
r  o  3 

o 


Eque  .ions  (87)  through  (89)  can  be  solved  to  give  the  desired 
initial  c<-  iditions  in  the  form 


dr  2B 
o  o 


Amplitude  (thousands  of  n  mi)  Amplitude  (thousands  of  n  mi) 
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Lunar  longitude,  9m  (0),  (deg) 
(a)  Lunar  effect 


Solar  longitude,  0  (0),  (deg) 
(b)  Solar  effect 


Fig.  12 — Amplitude  of  the  orbital  frequency  comoonent  of  8r  due  to  the 
sun  and  the  moon  as  functions  of  the  initio*  values  of  0  and  Qm 
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(91) 


60'  -  2A.  -  B  (92) 

o  io 

With  these  initial  conditions,  Kqs.  (76)  and  (77)  can  be  re¬ 
duced  to 


wher*  the  terns  of  the  homogeneous  solution  present  in  Eqe.  (76)  and 
(77)  are  eliminated  with  the  exception  of  a  small  bias  term  In  6r. 


-40- 


It  must  be  emphasized,  however,  that  these  linearized  solutions 
are  used  only  as  a  basis  for  determining  initial  conditions  and  not 
as  a  substitute  for  the  numerical  integration  of  the  actual  nonlinear 
equations . 

Equations  (90)  through  (92)  express  the  desired  initial  condi¬ 
tions  in  terms  of  A^,  A^,  and  Bq,  which  are  all  functions  of  the  ini¬ 
tial  values  of  0,  0  ,  and  0.  An  examination  of  Eqs.  (79)  through  (81) 
shows  that  A,,  A„,  and  B  are  each  summations  of  about  480  terns: 

XL  O 

these  expresions  can  be  programmed  for  evaluation  by  JOSS  so  that 

6r  ,  6r' ,  and  £6'  can  be  determined. 

O  0  o 

As  a  check  on  the  foregoing  analysis  the  following  case  is 
considered: 

Reference  Orbit 

Orbital  radius,  rQ  -  303,000  n  mi 

Orbital  period,  Tq  ■  48.53  days 

Orbital  velocity,  V  •  2761  ft/sec 

o 

Initial  Geometry 

Orbital  angle,  8(0)  •  0  deg 

Lunar  longitude,  6  (0)  -  40  deg 
Solar  longitude,  9(0)  -  0  deg 

Table  2  presents  a  comparison  of  the  reference  orbital  injec¬ 
tion  conditions  and  those  computed  from  Eqs.  (90)  through  (92)  as 
well  as  the  values  if  6r  and  3r.c,  which  should  result  in  each  case. 

OSC  Jm 

Table  2 

ORBITAL  INJECTION  CONDITIONS 
(Case  #1) 


Parameter  Reference  Computed 

Radial  increment,  4r  (n  mi)  0  21 
Horizontal  veloctiy  (ft/sec)  2761  2860 
Vertical  velocity  (ft/sec)  0  -61 


Oscillatory  amplitude,  6r09c  (n  mi)  22,849  i  0 
Steady-state  bias,  £rss  <n  mi)  -20,359  j  1522 


variation  (thousands 


Figure  13*  shows  the  resulting  orbital  variation  when  the  refer¬ 
ence  orbital  injection  conditions  are  used.  It  Is  seen  that  during 
the  Initial  part  of  the  orbit  the  bias  and  oscillatory  amplitude  are 

in  good  agreement  with  the  values  tabulated  for  ir„.  and  6r  ;  how- 

bb  osc 

ever,  as  time  goes  on,  the  amplitude  Increases  and  exceeds  the  50,000 
n  ml  stability  limit  after  187  days.  In  comparison.  Fig.  13b  shows 
the  corresponding  orbital  motion  when  the  computed  Injection  condi¬ 
tions  are  used.  The  resulting  magnitude  of  5r  remains  less  than 
17,000  n  ml  for  250  days  with  no  Indication  of  an  orbital  frequency 
oscillation.  From  Table  2  it  is  seen  that  this  marked  Improvement 


Init.  cend. 

Horiz.  vel. 

2761  ft/ t 

Rad.  Tel. 

C  ft/l 

Rad.  diip. 

0  n  mi 

Time  (days) 

(a)  Uncorrected 


Init.  eoad. 


Time  (days) 
(b)  Corrected 


rQ  =303,000  n  mi 
0m(O)=  40' 

,*  0° 


Fig,  12 — Improvement  in  orbital  stability  by  alteration 
of  initial  conditions 
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ls  achieved  for  on  increase  of  about  100  ft/aec  in  injection  velocity 
and  21  n  mi  in  injection  altitude. 

In  the  application  of  this  method  of  determining  injection  condi¬ 
tions,  there  are  certain  restrictions  on  the  selection  of  the  refer¬ 
ence  orbital  radius.  An  examination  of  Eqs.  (78)  through  (85)  shows 
that  values  of  or  nearly  equal  to  either  0  or  ±1  can  result 
in  very  large  terms  in  the  susmatlons.  Thus  it  is  necessary  that  r^ 
be  selected  so  that  the  corresponding  orbital  angular  rate  0q  is  in¬ 
commensurate  with  0  and  0  .  This  selection  is  discussed  in  more 

mo  o 

detail  in  Appendix  8. 


Orbital  Correction  Procedure 

It  must  be  recognized  that  the  initialization  method  described 
above  is  based  on  a  linearization  of  the  equations  of  motion  and  its 
usefulness  in  establishing  a  stable  orbit  depends  on  how  rapidly  the 
nonlinear  effects  begin  to  appear.  It  is  found  that  the  development 
of  a  divergent  oscillation  at  orbital  frequency  does  not  ordinsrily 
occur  until  at  least  a  year  after  orbital  injection,  and  stability 
in  the  sense  of  not  exceeding  50,00G  n  ml  for  6r  may  last  for  2  to 
3  years.  However,  it  appears  that  even  with  the  initialization  proce¬ 
dure,  instability  can  occur  in  less  than  five  years  unless  further 
corrective  action  is  taken. 

Orbital  correction  can  be  settle*  »d  by  a  modification  of  the 
method  used  to  determine  the  orbital  injection  condition  at  time  t 
equal  to  zero.  Zf  t^  Is  the  time  at  which  the  correction  is  to  be 
made,  then  the  linearized  solution  in  the  vicinity  of  t^  is  given  by 
the  expressions 
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sin  (6o  -  8X) 


+  driven  solutions 


(96) 


where 

6i  ■  Vi  (,7) 

while  Sr^t  dr^,  and  68^  are  the  variations  existing  at  tine  t^,  and 
the  values  of  Ao,  Ap  A^,  and  Bq  are  determined  by  Eqs.  (78)  through 
'81),  also  evaluated  at  tit*',  t. .  The  driven  solutions  are  the  same 
*.  those  appearing  In  Eqs.  (76)  and  (77),  since  they  are  Independent 
of  the  value  of  t^. 

Since  orbital  Instability  Is  characterised  by  a  buildup  of  an 

oscillation  in.  both  Cr  and  d«8/d6^  at  orbital  frequency,  the  current 

o 

amplitude  of  this  oscillation  can  be  used  as  a  criterion  for  initiat¬ 
ing  an  oxbltal  correction.  This  amplitude  at  time,  t^,  is  given  by 


When  this  amplitude  exceeds  some  preassigned  limit  the  necessity  of  an 
oxbltal  correction  la  indicated. 

Application  of  Eqs.  (90)  through  (92)  at  time  gives  the  fol¬ 
lowing  set  of  desired  initial  conditions: 


r 


o 


(99) 
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t-t. 


(100) 


59i  -  <2A1  -  Vt-t,  (10] 

As  before,  If  these  values  are  substituted  for  4r^/rQ,  dr^,/ro, 
and  50^  in  Eqs.  (95)  and  (96)  the  amplitude  of  the  orbital  frequency 
oscillations  and  the  steady-state  term  in  Eq.  (96)  vanish. 

Unfortunately,  this  procedure  is  only  possible  when  the  current 
value  of  6r^  happens  to  be  equal  to  the  computed  vnlue  of  dr^,  since 
step  changes  in  6r^  are  not  possible.  However,  an  impulsive  velocity 
change  can  produce  step  changes  in  both  dr^  and  50^.  By  changing 
these  two  quantities,  it  is  possible  to  reduce  the  coefficient  of 
sin  (0  -  9.)  in  Eqs.  (95)  and  (96)  to  zero  as  well  as  the  constant 

O  a 

term  in  Eq.  (95) .  The  corrected  values  of  5r^  and  69^  are  given  by 


*“A2*t-t, 


(102) 


(103) 


Substitution  of  Eqs.  (99)  through  (103)  into  Eqs.  (95)  and  (96) 
results  in  the  following  form  of  the  linearized  solution  in  the  vicin¬ 
ity  of  tj^: 


-i  _  cos  (9o  -  ej}  +  driven  solutions 


(104) 
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Examinatioo  of  Eq*.  (104)  and  (105)  shows  that  if  the  current 

value  dr^  equals  the  desired  value  drD  computed  front  Eq.  (99),  then 

the  linearised  solutions  are  similar  to  those  given  previously  in 

Eq».  (93)  and  (94)  with  no  orbital  frequency  oscillation  and  no  con-* 

st ant  term  in  ddd/dfl  .  If  dr,  and  dr_  are  not  equal,  the  resulting 

o  in 

value  of  dr  Is  equal  to  thiir  difference, 
osc 

The  horizontal  and  radial  velocity  changes  required  to  accom¬ 
plish  this  correction  lire  given  by  the  relations 


*VR“ 


(107) 


The  correction  time  t^  can  be  determined  by  computing  the  val¬ 
ues  of  dr^  and  drp  as  functions  of  time  for  several  days  in  advance, 
thus  permitting  the  selection  of  a  value  of  t^  such  that  dr^  and  dr^ 
are  equal. 

This  procedure  is  applied  using  values  of  303,444  n  mi  for  r  , 

v 

158. 8  deg  for  8^(0),  0  deg  for  0(0),  and  10,090  n  mi  for  the  threshold 
value  of  dr  .  The  resulting  variation  of  dr  as  a  function  of  time 
is  shown  in  Fig.  14,  where  the  positions  at  which  corrections  ere  ap¬ 
plied  ere  indicated  by  a  vertical  line.  After  each  correction  the 
solid  curve  represents  the  corrected  notion,  while  the  dotted  curve 
shows  the  motion  thet  results  if  no  correction  is  made. 

The  initial  conditions  at  orbital  injection,  as  well  as  tlv  e 
for  the  reference  orbit,  are  shown  In  Table  3,  while  Table  4  shows 
the  details  of  the  four  corrections  applied  including  the  tine  of 
application,  the  horisootel  and  vertical  components  of  the  velocity 
increment,  its  total  magnitude  and  the  accuaulated  magnitude  of  all 
of  the  increments. 

Prom  Fig.  14  and  Table  4  it  is  seen  thet  for  s  total  velocity 
correction  of  342.6  ft/eec,  the  magnitude  c *  dr  remains  less  then 


Corr 


ig  orbital 
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Orb  leal  radius  (n  mi) 
Horizontal  velocity  (ft/aec) 
Vertical  velocity  (ft/3ec) 


Computed 

Reference 

305,762.4 

2790.09 

9.80 

303,444.0 

2758.59 

0 

Table  4 

ORBITAL  CORRECTIONS 
(Case  #2) 


Velocity  Increment 


Time 

Number  (days) 


235.5 

501.5 
909.0 

1413.5 


Horizontal 

(ft/aoc) 


Vertical 

(ft/aec) 


Total 

(ft/aec) 


65,3 

117.3 

69.9 

90.1 


Accumulated 

(ft/aee) 


65.3 

182.6 

252.5 

342.6 


25,000  n  ml  for  a  period  of  five  years.  This  represents  an  expendi¬ 
ture  of  68.5  ft/aec  ptr  year. 

The  choice  of  a  threshold  value  of  10,000  n  nd  for  the  value  of 
dr  at  which  orbital  correction  1&  initiated  may  seem  low  in  view 

OS  C 

of  the  50,000-n  mi  stability  limit  on  fr  specified  earlier.  However, 
an  examination  of  Table  4  shows  that  the  major  part  of  the  correction 
is  Involved  in  reducing  the  magnitude  of  the  radial  velocity.  If  the 
oscillation  at  orbital  frequency  reaches  an  amplitude  of  50,000  n  mi, 
the  maximum  radial  velocity  is  of  the  order  of  450  ft/sec  and  occurs 
when  6r  is  near  zero.  Since  the  correction  also  occurs  when  6r  is 
small,  a  considerable  fraction  of  this  450  ft/sec  is  required  as  a 
corrective  velocity  impulse  to  reduce  this  radial  velocity.  This  is 
more  than  the  accumulated  value  of  342.6  ft/sec  shown  in  Table  4  for 
the  four  corrections  applied  in  that  case.  This  magnitude  of  correc¬ 
tive  impulse  in  itself  is  not  important,  if  the  buildup  of  the 
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oscillation  Is  linear.  Under  such  conditions,  it  is  immaterial  whether 
the  orbital  correction  is  achieved  by  small  velocity  inpulses  at  short 
Intervals  or  large  impulses  at  longer  intervals,  since  the  accumulated 
impulse  is  the  same.  However,  It  is  found  that  the  buildup  is  not 
linear,  and  thus  it  is  more  economical  to  correct  more  frequently.  In 
addition,  this  procedure  results  in  maximum  values  for  dr  considerably 
less  than  the  50,000-n  mi  limit  originally  specified. 
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IV.  CONCLUSIONS 


On  the  basis  of  the  foregoing  results  and  discussion,  the  follow¬ 
ing  specific  conclusions  can  be  stated  regarding  the  motion  of  earth 
•7  Jte1  '  ‘  it  trans lunar  distances  from  the  earth. 

o  The  existence  of  stable  earth  orbits  in  translunar  space  seems 
more  likely  at  distances  between  300,000  \>  ^00,000  n  mi  from 

the  earth,  since  in  this  range  the  maximum  gravitational  dis¬ 
turbance  due  to  the  sun  and  the  moon  is  minimized  relative  to 
the  earth *s  attraction. 

o  A  value  of  orbital  radius  closer  to  the  300,000-n  mi  limit  ap¬ 
pears  to  be  preferable,  since  the  average  gravitational  dis¬ 
turbance  due  to  the  moon  is  considerably  less  than  its  maximum 
during  a  clcse  approach  of  the  satellite  to  the  moon. 

o  The  steady-state  motion  of  the  orbital  plane  of  &  translunar 
earth  satellite  is  a  slow  regression  at  constant  inclination 
to  the  plane  of  the  ecliptic  with  small  oscillatory  variations 
in  regression  rate  and  inclination  superposed.  Since  the  or¬ 
bital  planes  of  near-earth  satellites  are  subject  to  similar 
motions,  this  can  hardly  be  regarded  as  orbital  instability. 

o  The  period  of  orbital  regression  increases  nonotonically  with 
inclination  to  the  ecliptic,  ranging  from  five  years  at  0  deg 
inclination  to  infinity  at  90  deg  and  being  relatively  insen¬ 
sitive  to  the  value  of  orbital  radius  in  the  300,000-  to 
500,000-n  ml  range. 

o  The  orbital  inclination  to  the  equate rial  plane  is  not  neces¬ 
sarily  constant  and  nay  vary  aa  much  as  45  deg  in  intervale 
as  short  as  three  or  more  years. 

c  Tho  only  conditions  under  which  the  orbital  inclination  to  the 
equator  can  remain  constant  are  those  for  an  orbital  plane  in 
the  ecliptic  or  perpendicular  to  It. 

o  An  orbit  in  the  ecliptic  maintains  a  constant  equatorial  in¬ 
clination  of  23*27' ,  while  the  corresponding  constant  inclina¬ 
tion  for  an  orbit  perpendicular  to  the  ecliptic  may  have  any 
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value  between  66*33'  and  90  deg  depending  on  the  fixed  orien¬ 
tation  of  Its  line  of  nodes  in  the  ecliptic, 
o  The  choice  of  an  orbital  plane  perpendicular  to  the  ecliptic 
tends  to  reduce  the  magnitude  of  the  disturbance  due  to  the 
noon  both  by  decreasing  the  probability  of  a  near  approach  to 
the  noon  and  by  reducing  its  duration  if  it  does  occur.  This 
should  in  turn  improve  the  characteristics  of  the  motion  of 
the  satellite  in  the  orbital  plane, 
o  Orbital  instability  can  develop  in  connection  with  the  in¬ 
plane  motion  of  the  satellite  and  is  characterized  by  a  di¬ 
vergent  orbital  frequency  oscillation  in  the  radial  displace¬ 
ment  from  a  nominal  circular  orbit, 
o  By  a  cut  and  try  approach,  initial  conditions  can  sometimes 
be  found  so  that  the  amplitude  of  the  radial  oscillations  re¬ 
mains  less  than  50,000  n  ml  for  periods  as  long  as  five  years, 
o  The  orbital  Injection  conditions  can  be  determined  in  terms 
of  the  initial  positions  of  the  satellite,  the  sun,  sad  the 
moon  relative  to  the  earth  so  that  the  development  of  the  os¬ 
cillatory  instability  is  delayed  for  two  to  three  years, 
o  It  is  also  possible,  if  instability  begins  to  develop,  to  de¬ 
termine  the  magnitude  and  time  or  application  of  an  impulsive 
velocity  correction  to  reduce  the  oscillation, 
o  Application  of  this  method  of  orbital  correction  indicates 
that  it  is  possible  to  keep  the  amplitude  of  the  radial  os¬ 
cillation  below  25,000  n  *<L  for  a  period  of  five  years  at  an 
expenditure  of  68.5  ft/sec  per  year  in  velocity  impulse. 

To  susmarlze  these  conclusions.  It  appears  that  although  orbits 
can  be  determined  that  have  satisfactory  stability  over  a  five-ye*r 
interval,  the  actual  realization  of  such  an  crbit  may  be  highly  sen¬ 
sitive  to  the  accuracy  with  which  the  orbital  injection  conditions  are 
achieved.  As  a  result,  it  seems  advisable  to  provide  an  orbital  cor¬ 
rection  capability  to  insure  the  desired  performance. 


* 
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Append  lx  A 

SERIES  EXPANSION  OF  THE  SOLAR  AND  LUNAR 
DISTURBING  ACCELERATIONS 


EQUATIONS  OF  MOTION 


A  linearized  form  of  the  equations  for  the  in-plane  motion  of  a 
satellite  with  a  90-dejj  orbital  inclination  to  the  ecliptic  is  devel¬ 
oped  in  the  body  of  this  report  as  Eqs.  (66)  and  (67)  in  the  font 


o 


(A-l) 


(A-2) 


where  aSr,  aff,  ,  afi9  are  the  normalized  radial  and  tangential  dis¬ 
turbing  accelerations  due  to  the  sun  and  the  moon.  These  quantities 
are  expressed  as 


*5r 


2 

cos 


e) 


(A-3) 
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For  use  in  Eqo.  (A-l)  and  (A-2)  it  is  more  convenient  to  express  these 
normalized  accelerations  as  a  summation  of  sinusoids  in  terms  of  the 
Independent  variable  8q.  This  is  done  in  the  following  subsections. 


SOLAR  ACCELERATION 

The  desired  expansion  of  the  solar  terms  is  as  follows: 


e;  r  i  i 
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which  can  be  written  as 
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where  the  nonzero  values  of  are 
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which  be cooes 
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where 
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Equation  (A-l4.)  can  be  expressed  as  a  binomial  expansion  in  the  font 
of  an  infinite  series  as 


(M3 


Ka  CM®  6  cos°  * 


(A-16) 


which  can  be  approximated  by  truncating  the  series  after  the  first  H 

terns  as 
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cos  b  cos  0 


(A-17) 


Since  the  -value  of  K  for  a  300,000-n  mi  orbit  is  0.4356,  the  con¬ 
vergence  of  the  series  Is  not  rapid,  particularly  if  6  and  are  both 
equal  to  either  0  or  180  deg.  As  a  result,  a  value  of  at  lease  100 
for  N  is  required  to  give  a  satisfactory  approxlaatlon  to  the  actual 
function. 

Equation  (A-17)  can  be  further  expanded  by  aaking  use 

sj  H 

Fourier  expansion  of  cos  0  and  cos  0  in  the  fora 
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These  coefficient*  are  related  by  the  following  recursion  forwcl**: 
(*  If  .-0 
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Equation  (A-17)  can  now  become 
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The  product  of  the  two  inner  summations  can  be  expressed  In  terms 
of  sums  and  differences  of  thv*.  angles  10  and  j0ffl  so  that  Eq.  (A-24) 


becomes 
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By  revTsrstn,*:  the  order  of  the  summations  In  Eq.  (A-25)  it  can  be 
,'.ewr<  .n  as 
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where 


3/2  N 


■©  s* 


2  (nt) 


21  KnB 

2  K  nij 


(A-28) 


and  n  is  the  highest  harmonic  of  the  6  and  6  motion  necessary  to 
o  m  ^ 

give  an  accurate  expansion  of  the  function  (a  /v  )  .  It  is  found 
*  mm 

that  a  value  of  20  for  nQ  is  satisfactory e 

Equation  (A-27)  is  tne  desired  expansion  for  the  first  term  of 

a  „  as  given  In  Eq.  (A-5) ,  and  can  he  used  to  obtain  a  similar  expan- 
mr 

8 ion  of  the  second  term  as  follows t 
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Substitution  of  Eqs .  (A-25)  and  (A-30)  in  Eq.  (A-5)  gives 
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Equation  (A-32)  is  the  desired  expansion  of  amr  as  a  summation 
of  sinusoidal  terms  with  arguments  that  are  st-wa  and  differences  of  the 
first  n  harmonics  of  6  and  8^. 

The  quantity  amQ  can  be  expanded  in  a  simitar  manner,  its  first 
term  from  Eq.  (A-6)  being  given  by 
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ACCURACY 

In  the  development  of  the  series  expansion  for  a  and  a  a  above, 

mr  mo 

the  binomial  expansion  in  Eq,  (A-17)  is  terminated  after  100  terms. 

In  addition,  harmonics  cf  0  ar.d  0  above  the  twenty-first  are  neglected 

m 

in  the  Fourier  series  for  each  of  the  terms  of  the  binomial  expansion. 

As  a  check  on  the  accuracy  of  this  truncation,  Table  5  gives  a  compari¬ 
son  between  the  values  of  a  „  computed  from  the  series  of  Eq,  (A-32) 
and  from  the  exact  expression  given  by  Eq.  (A-5) .  Table  6  represents 
a  similar  comparison  between  Eqs.  (A-6)  and  (A-37)  for  amg.  These  eval¬ 
uations  are  made  for  a  300,000-n  mi  orbit  in  the  vicinity  of  a  close 

approach  where  0  and  0  are  both  zero,  and  even  under  this  extreme 
o  m 

condition,  the  series  expansion  is  in  error  by  less  than  0.5  percent. 
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Table  5 
ACCURACY  OF  a 

mr 


9  (deg) 

i 

0ffl  (deg) 

a 

t 

nr 

Series 

Exact 

0 

0 

-0.154226 

-0.154966 

1 

1.75 

-0.152744 

-0.153405 

2 

3.50 

-0.148458 

-0.148919 

3 

5.25 

-0.141819 

-0.142047 

4 

7.00 

-0.133474 

-0.133523 

5 

8.75 

-0.124143 

-0.124107 

6 

10.50 

-0.114499 

-0.114461 

7 

12.25 

-0.105080 

-0.105080 

8 

14.00 

-0.096248 

-0.096284 

9 

15.75 

-0.088197 

-0.088242 

10 

17.50 

-0.080983 

-0.081015 

15 

26.25 

-0.055594 

-0.055585 

20 

35.00 

-0.041447 

-0.041450 

25 

43.75 

-0.032389 

-0.032388 

30 

52.50 

-0.025632 

-0.025631 

Table  6 


ACCURACY  OF  a  . 

mO 


6  (deg) 

8ffl  (deg) 

a 

a 

i6 

Series 

Exact 

0 

0.00 

0.000000 

-0.000000 

1 

1.75 

-0.004509 

-0.004535 

2 

3.50 

-0.008634 

-0.008670 

3 

5.25 

-0.012064 

-0.012091 

4 

7.00 

-0.014609 

-0.014617 

5 

8.75 

-0.016217 

-0.016211 

6 

10.50 

-0.016955 

-0.016947 

7 

12.25 

-0.016965 

-0.016966 

8 

14.00 

-0.016420 

-0.016429 

9 

15.75 

-0.015483 

-0.015497 

10 

17.50 

-0.014292 

-0.014303 

15 

26.25 

-0.007352 

-0.007349 

20 

35.00 

-0.001749 

-0.001750 

25 

43,75 

+0.001865 

+0.001865 

30 

52.50 

+0.003839 

+0.003839 

%  f  ■^Twrf  I'nrs^v.'y^r  'j- 
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AppendJ.x  B 

SELECTION  OF  A  REFERENCE  ORBITAL  RADIUS 


In  the  body  of  the  report,  a  method  is  described  for  the  determi¬ 
nation  of  the  orbital  injection  conditions  dr  ,  Sr' ,  and  60'  in  terms 

0  o  o 

of  three  summations  A^,  and  Bq  given  by  Eqs.  (79)  through  (81).  An 
examination  of  these  equations  shows  that  if  any  of  the  values  of 
or  are  near  zero  the  corresponding  term  in  Aq  tends  toward  infinity. 
Similarly,  if  or  approach  either  +1  or  -1,  the  corresponding 
terms  in  A^  and  A^  can  become  infinite.  In  either  case  the  injection 
conditions  computed  from  Eqs.  (90)  through  (92)  are  not  valid.  The 
reason  for  this  behavior  can  be  seen  by  examining  the  driven  solution 
for  a  single  frequency,  where 
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where 


16(0)  +  j6  (0) 
in 


(B— 3) 


Ordinarily  the  first  term  in  Eq.  (B-l)  is  included  in  the  summation  for 

A  ,  while  the  coefficients  of  cos  6  and  sin  '  are  included  in  A.  and 
o  o  o  1 

respectively.  Likewise,  the  first  term  of  Eq.  (B-2)  is  part  of  the 
suamatlon  for  Bq.  However,  sa  approaches  0,  the  above  solution  be¬ 
comes  indeterminate  and  Eqs.  (B-l)  and  (B-2)  approach  the  forms 


T-  ■  ’uV*  “*  *U 


(B-4) 
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Similarly,  if  approaches  either  +1  or  -1,  the  solutions  become 
«r  _  (p11  *  2V  ...  . 


2 - 8  sin  (0  ±  #  ) 


(B-6) 


IT -  -(P^  *  2qlj)0  sin  (6  ¥ 
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Since  the  values  of  p^  and  q^  are  much  less  than  1,  the  resulting  so¬ 
lutions  for  6r/rQ  *ad  dd6/d9Q  in  Eqs.  (B-4)  through  (B-7)  diverge  from 
0  very  slowly  even  though  in  each  case  the  solutions  are  the  difference 
of  two  large  quantities.  As  a  result,  it  does  not  seem  reasonable  to 
include  one  of  these  large  terms  in  AQ,  A^,  A^,  or  BQ  to  be  cancelled 
by  the  choice  of  Initial  conditions  when  this  leaves  a  nearly  equal 
term  uncompensated  in  the  driven  solution. 

However,  these  resonance  conditions  can  be  avoided  by  a  suitable 
choice  of  the  reference  orbital  radius.  From  Eqs.  (24),  (74)  and  (75) 
the  dependence  of  and  on  Tq  is  given  as 
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For  the  range  of  values  of  rQ  considered  in  this  report  and  the 
values  of  i  and  j  involved,  never  becomes  equal  to  either  0  or  ±  1. 

This  is  not  the  case  with  w.  .  as  can  be  seen  from  Fig.  15  where 
the  values  of  i  are  plotted  horizontally  and  those  of  -j  vertically. 

The  points  on  the  i,j  plane  represent  the  combination  of  integral  val¬ 
ues  of  i  and  j  for  which  an  exists.  The  three  parallel  lines  origi¬ 
nating  at  -1,  0,  and  4-1  on  the  1  axis  represent  the  contours  in  the  i,j 
plane  along  which  is  equal  to  -1,  0,  and  +1,  respectively,  and  are 
given  by  Eq.  (B-8)  as 
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i  *  -  l 
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e 

,  _  .mo 
i  -  -ij- 

o 

i  -  ~^r  ♦ 1 

o 

where  the  slope  of  these  lines  is  a  function  of  through  Eq.  (B-10). 

To  avoid  the  resonance  conditions  it  is  necessary  to  select  a  value  of 
rQ  so  that  none  of  the  three  parallel  lines  past  through  any  of  the  fre¬ 
quency  points  in  the  i,j  plane.  A  simplified  version  of  this  determina¬ 
tion  is  shown  in  Fig.  16  where  all  of  the  integer  combinations  of  i  and 
j  are  indicated  as  frequency  points  instead  of  just  those  for  which  i 
and  j  are  both  even  or  both  odd.  Thus,  if  a  straight  line  from  the 
origin  misses  all  of  these  frequencies,  then  Eqs.  (8-11)  through  B-13) 
are  satisfied  and  the  intersection  of  this  line  with  the  vertical  r 

o 

scale  determines  an  acceptable  reference  orbital  vadius. 


(3-11) 

(B-12) 

(B-13) 
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Pig.  16 — Determination  of  reference  radius  to  avoid  resonance 


Orbital  radiui  (thousands  of  n  mi) 
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At  first  glance  it  may  appear  that  it  is  impossible  to  avoid  coming 

close  to  some  of  these  frequency  points.  However,  the  values  of  p_^  and 

q^j,  which  occur  in  the  umerators  of  the  terms  in  Eqs,  (B-l)  and  (B-2) , 

are  much  less  than  unity  and  the  resulting  resonance  is  very  sharp.  As 

an  example,  the  value  300,000  n  mi  for  rQ  turns  out  to  be  a  poor  choice, 

since  the  corresponding  ratio  of  9  /6  is  equal  to  1.74998782,  which 

mo  o 

is  a  good  approximation  to  7/4.  An  examination  of  Fig.  16  shows  that  a 
line  from  the  origin  to  the  300,000-n  mi  position  on  the  orbital  radius 
scale  passes  very  close  to  the  points  (i  ■  7,  j  *  -4),  (i  ■  14,  j  ■  -8), 
and  (i  »  21,  j  *  -12) .  This  indicates  that  the  following  resonances  are 
possible 

“8,  -4  "  1 


U),  ,  *  “1 

6,  -4 


“14.  -8  ’  0 


“20,  -12  ■  -1 


In  Table  7.  the  amplitudes  and  c^  of  the  last  terms  in  Eqs. 
(B-l)  and  (B-2)  are  tabulated  for  all  terms  with  a  value  of  l^jl  lees 
than  lo5  for  an  orbital  radius  of  300,000  n  mi.  It  is  seen  that  all 
four  of  the  predicted  resonances  t,  but  the  resulting  amplitudes  are 
much  smaller  for  the.  higher  values  of  i  and  j.  This  is  primarily  due 
to  the  fact  that  the  quantities  p^  and  which  appear  in  the  numer¬ 
ators  of  a^  and  c^,  descrease  rapidly  as  i  and  j  become  larger.  In 
comparison,  Table  8  gives  the  corresponding  amplitudes  when  the  value 
of  rQ  is  shifted  from  300,000  to  303,000  n  mi  and  shows  no  significant 
resonance. 

Thus,  by  a  relatively  small  shift  in  the  reference  orbital  radius, 
the  resonances  can  be  avoided  and  the  method  described  for  computing 
orbital  injection  conditions  can  be  used. 


x 
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Table  7 


RESONANCE  EFFECTS 
(rQ  *  300,000  n  mi) 


i 

j 

Kji 

Cij 

1 

-1 

0.74999 

0.00612 

0.00052 

2 

-2 

1.49998 

0.00594 

0.01358 

3 

-1 

1.25001 

0.00149 

0.00507 

4 

-2 

0.50002 

0.00583 

0.00847 

5 

-3 

0.24996 

0.00934 

0.01490 

6 

-4 

0.99995 

18.68356 

37.36655 

8 

-4 

1.00005 

2.62238 

5.24508 

9 

-5 

0.25006 

0.00135 

0.00195 

10 

-6 

0.49993 

0.00077 

0.00131 

11 

-7 

1.24991 

0.00033 

0.00071 

13 

-7 

0.750C9 

0.00013 

0.00020 

14 

-8 

0.00010 

0.44030 

0.66043 

15 

-9 

0.74989 

0.00011 

0.00021 

16 

-10 

1.49988 

0.00002 

0.00004 

17 

-9 

1.25011 

0.00001 

0.00002 

18 

0.50012 

0.0U002 

0.000C2 

19 

0.24987 

0.00002 

0.00004 

20 

IS 

0.99985 

0.03023 

0.06046 

Table  8 


RESONANCE  EFFECTS 

(r  *  303,000  n  mi) 
o 


i 

J 

Kjl 

, 

“ij 

clj 

1 

-1 

0.77630 

0.00617 

0.00065 

3 

-1 

1.22370 

0.00181 

0.00569 

4 

-2 

0.44739 

0.00615 

0.00888 

5 

-3 

0.32891 

0.00738 

0.01204 

6 

-4 

1.10521 

0.00725 

0.01497 

8 

-4 

0.89479 

0.00151 

0.00271 

9 

-5 

0.11848 

0.00268 

0.00394 

10 

-6 

0.65782 

0.00074 

0.00133 

11 

-7 

1.43412 

0.00014 

0.00033 

12 

-6 

1.34218 

0.0PQ02 

0.00008 

13 

-7 

0.56588 

0.00011 

0.00016 

14 

-8 

0.21043 

0.00021 

0.00032 

15 

0.9S673 

0.00131 

0.00261 

17 

1.01327 

0.00019 

0.00038 

18 

IMIk 

0.23897 

0.00003 

0. 00004 

19 

0.53933 

0.00001 

0.00002 

20 

-12 

1.31564 

o.ooooi 

0.00002 

21 

-11 

1.46067 

0.00000 

0.00001 
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Appendix  C 

JOSS*  COMPUTING  PROGRAMS 


This  appendix  includes  the  programs  used  for  the  Runge-Kutta  Inte¬ 
gration  of  the  Orbital  Equations  and  for  the  computation  of  Orbital 
Injection  and  Correction.  Following;  each  program  is  a  sample  of  the 
output,  an  Identification  of  JOSS  terms,  and  a  description  of  the  in¬ 
dividual  parts  of  the  urogram. 


JOSS  is  the  trademark  and  service  mark  of  The  Rand  Corporation 
for  its  computer  program  and  services  using  that  program. 


v 


-71- 


Runge-Kutta  Integration  of  Orbital  Equations 

*****************;<**************************** 


1.001  ^Initiation  (Do  part  1.) 

1.03  Demand  A(0)  ea  "Reference  Orbital  Radius  (n  mi)". 

1.04  Demand  A(l)  as  "Starting  Time  (days)". 

1.05  Demand  A(2)  as  "Terminal  Time  (days)". 

1.06  Demand  A(3)  as  "Orbital  Inclination  to  Ecliptic  (deg)". 

1.07  Demand  A(4)  as  "Computing  Step  Size  (days)". 

1.08  Demand  A(5)  as  "Printout  Step  Size  (days)". 

1.09  Demand  A(61)  as  "Initial  Orbital  Central  Angle  (deg)". 

1.10  Set  A(62)sA(51)*arg(-l,0)/180. 

1.11  Demand  A(6)  as  "Initial  Lunar  Longitude  (deg)". 

1.12  Demand  A(7)  as  "Initial  Solar  Longitude  (deg)". 

1.13  Demand  A(52)  as  "Initial  Longitude  of  Lunar  Perigee  (deg)". 

1.14  Demand  z(l,0)  as  "dr/R". 

1.15  Demand  z(ltl)  as  "d(dr/R)/dQ". 

1.16  Demand  z(2,0)  as  "dq  (rad)"  if  A(1)*0. 

1.17  Set  z(2,0)*0  if  A(1)«0. 

1.18  Demand  z(2,l)  as  "d(dq)/dQ". 

1.19  Demand  z(3,0)  as  "dq(M)  (rad)"  if  A(1)*0. 

1.20  Set  z(3#0)=G  if  A(1)=0. 

1.21  Demand  z(4t0)  as  "dq(S)  (rad)"  if  A(1)*0. 

1.22  Set  z(4,0)=0  if  A(1)*0. 

1.23  Do  part  2. 

1.24  Do  part  3. 

1.25  Do  part  5. 

2.001  *Evaluation  of  Parameters 
2.01  Set  A(8)*82.27. 

2.02  Set  A(9)*arg(-1,0). 

2.03  Set  A(10)*A(9)/180. 

2.04  Set  A(ll)«238857*5280/6080. 

2.05  Set  A( 12 )*9. 3» 10*7*5280/6080 • 

2.06  Set  A( 13)*. 0549. 

2.07  Set  A(14)». 01674. 

2.08  Set  A(15)*»0019437. 

2.09  Set  A( 16 )*[ 281+33/60+53/3600] » A( 10) . 

2.10  Set  A(17)»A(0)*6080, 

2.11  Set  A(18)*.22998*sqrt[(A(8)-l)/A(8)»A(ll)*3/A(0)*3]. 

2.12  Set  A(19)*A(18)/86400. 

2.13  Set  A( 20 )*A( 17 ) « A( 19 ) « 

2.14  Set  A(21)«A(17).A(19}*2. 

2.15  Set  A(22)*2*A(9)/A(18). 

2.16  Set  A(23)*. 22998. 

2.17  Set  A(24)«.0172. 

2.18  Set  A(25)«(A(23)/A(18))*2/Av8). 

2.19  Set  A(26)*(A(24)/A(18))*2. 

2.20  Set  t*A( 1 ) • 

2.21  Set  q*A(l)*A(18). 

2.22  Set  A(30)*0. 

2.23  Set  A( 31)*A(23)/A(18). 

2.24  Set  A(32)*A(24)/A(18). 
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2.25  Set  A(33)«A(15)/A<18), 

2.26  Set  A(34)*A(3i)-A(33). 

2.27  Set  A(53)*A(52)»A(10). 

3.001  *Fozaat 
3.01  Page. 

3.02  Type  A(0)  in  fora  1. 

3.03  Type  A(3)  in  fora  2. 

3.04  Type  A(18)  in  for*  3. 

3.05  Type  A(22)  in  fora  4. 

3.06  Type  A(21)  in  fora  S. 

3.07  Type  A(4)  in  fora  6. 

3.08  Type  [A(61)<0:A(61)+360;A(61)3  in  fora  18. 

3.09  Type  [A(6)<0:A(6)+360$A(6)3  in  fora  22. 

3.10  Type  [A(7)<0sA(7)+360}A<7)3  in  fora  23. 

3.11  Type  A(52)  in  fora  24. 

3.12  Line. 

3.13  Do  part  4  if  t*A(l). 

3.14  Type  fora  30,  foza  7, fora  8, form  31. 

4.001  ^Initial  Condition  Printout 
4.01  Type  t  in  fora  14. 

4.02  Line. 

4.03  type  x(l  ,0 )  ,z(  1 ,0 } *A(0 )  in  fora  10. 

4.04  Type  x(l,l),x(l,ljVA(20)  in  fora  11. 

4.05  Type  x(2,0),z(2,0)/A(10)  in  foza  12. 

4.06  Type  z(2,l),(z(2,l)+l)«r(x)oA(20)  in  fora  13. 

4.07  Type  *(3,0)  in  fora  20. 

4.08  Type  x(4,0)  in  foza  21. 

4.09  Line. 

4.10  Done  if  A(30)*l. 

4.11  Do  part  3  if  ft-A(2)3*agn(A(4))<0  and  t*A(l). 

5.001  ^Control  of  Coaputation  aui  Printout 
5.01  Do  part  6  if  f*>Ct/A<5)]»0. 

5.02  Set  A(29)*dlejC(t-A(2))*sgn(A(4))20,CQnj[A(30)Bl,fprt/A(5)3*033. 
5.03  type  foza  31,__  if  A(29)  or  $*43. 

5.04  Do  part  4  if  A(29)  or  $*45. 

5.05  Done  if  A(29). 

5.06  Do  part  7  for  i*9.01(  .01)9,03. 

5.07  Set  t»t+A<4). 

5.08  Set  q*qth. 

5.09  To  atep  5.01. 

6.001  ^Printout 

6.01  Set  i(l)*z(l,0)*A(0). 

6.02  Set  i(2)«z(2,0)/A(10). 

6.03  Set  i(3)*M(q,z)/A(10)+tv(lf(q,z)<G)*360. 

6.04  Set  i(4)*S(q,x)/A(10)+tv(S(q,z)<0)*360. 

6.05  Set  i(5)«arg[x(H,q,z),aqrt(l-x(M,q,z)*2)3/A(10). 

6.06  Set  i(6)*z(l,l)*A(20). 

6.07  Set  i(7)*b(q,z). 

6.08  Set  i(8)«r(zMN(x,q,z)-A(25)*d<q)*p(q,z)-A(26)«e{q)3. 

6,09  Set  i(9)*(z(2,l)tl)‘r(z)»A(20). 


f 
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6.10  Set  j(l)*D(q,z). 

6.11  Set  j(2)*N(y,q,z). 

6.12  Type  t,i(l)  ,i(6),i(7),i(8)  ,i(9),j(l),;3(2),i(2)  ,i(3),il4),i(5  j  in  fom  15. 

7.001  ^Variable  Selection 
7.01  Do  part  8  for  k*l(l)4. 

8.001  *Derivative  Selection 
8.01  Do  step  i  for  j*0,l  if  k<3. 

8.02  Do  step  1  for  j*0  if  k23. 

9.001  *Runge~Kutta  Integration 
9.01  Set  K(l,k,j)*h«J(q,z). 

9.02  Set  Z(k, j )*>z(k,  j )+K(l,k, j )/2. 

9o03  Set  K(2,k,j)=h»J(q+h/2,Z). 

9.04  Set  Z<k»j)«z(k,j)+K(2,k,j)/2. 

9.05  Set  K(3,k,j)*W(q+h/2,Z). 

9.06  Set  Z(k , j )«z(k , j )+K( 3 ,k , j ) . 

9.07  Set  K(4,k,j)«W(q+h,Z). 

9.08  Set  z(k,j)»z(kfj)+[K(i,k,j)+2«K<2,k,j)*2«K(3,k,j)tK(4,k,j)3/6. 


B(q,z) l 
C(q,z): 
D(q,z): 
F(q,z): 
J(q,z)i 
L(q.z): 
M(q,z) 
N(y,q.z) 
o(y,q,z) 

Q(q.*) 

R(q,z) 

S(q,z) 

T(q,z) 

a 

b(q,z) 

c(x) 

d<q) 

e(q) 

h 

p(q.z) 

r(z) 

*(x) 

x(»,q,z) 

y(«,q,z) 


[j*l:z(l,j+l);b(q,z)3 

[j*l:z(2,j+l);D(q,z)3 

-2»(l*z(2,l))/r(z)«z(l,l)+H(y,q,z) 

A(31)«LCl+A(13)*c(K(q,z)-A(33).q-A(53))]*2/(l-A(13)*2)*1.5-l] 

[k*l:B(q,z)sk*2:C(q,z);k*3:F(q,z)}k*4jL(q,z)3 

A(32).[[l+A'14).c(S(q,z)-A(16))]*2/(l-A(i4)*2)*1.5-l] 

fpC(A(31)»q+A(6)»A(10)+z(3,0))/2/A(9)]»2*A(9) 

A(25)*d(q)«Cp(q,z)-l3/R(q,z)*y(M,q,z)+0(y,q,z) 

3*A(26)*e(q)*x(S,q,z)*y(S,q,z) 

fp[(A(62)+q+z(2,0))/2/A(9)3*2»A(9) 

A(0)*r(z)/A(ll)*d(q)*(l/3) 

fpr(A(32)«q+A(7).A(10)+z(4,0))/2/A(9)3*2*A(9) 

2«z(2,l)+z(2ti)*2-A(25)»d(q)»p(q,z)+N(x,q,z)-A(26)*e(q) 

A(3)*A(10) 

z(l,0)/r(z)*2»[3+3*z(l,0)+z(l,0)*23+r(z)»T(q,z) 

coe(x) 

CCl+A(13)*c(M(q,z)-A(33)«q-A(53))3/Cl-A(13)*233*3 

CCl+A(i4).c(S(q,z)-A(16))3/[l-A(14)*233*3 

A(4)»A(18) 

[l+R(q,z)*2-2»R(q,z)»x(M,q,z)3*(-1.5) 

l+z(l,0) 

sin(x) 

c(Q(q,z))«c(K(q,z))+s(Q(q,z))*s(«(q,z))*c(a) 

»s(Q(q,z))»c(*(q,z))+c(Q(q,z))*s(B(q,z))*c(a) 


-74- 


Form  1: 

Orbital  Radius  n  ** 

Form  2: 

Orbital  Inclination  to  tha  Ecliptic  _ . _ dag* 

Form  3: 

Kaan  Orbital  Rata  *  rad/solar  day 

Form  4j 

Orbital  Period  *  days 

Form  5: 

Normal  Orbital  Acceleration  *  ft/sec/aec 

Form  6: 

Computing  Stop  Size  ,  days 

Form  7: 

|  t  |  dr  ! V(r)|  a(r)  |F<r)/m|V(q)|  a(q)  |F(q)/m|  dq  |Q(H)|Q(S)|Q(sm)| 


Form  8: 

I  days  I  n  mi  !ft/s|  I  I  ft/a  I  I  I  deg|  degl  dag  I  dag  I 

Form  10: 

dr/R  * . .  dr  * _ n  mi 

Form  11: 

d(dr/R)/dQ  * .  V(r)  * _ ft/a 

Form  12: 

dq  (rad)  *  •••••••«••••••  dq  *  _ * _ dag 

Form  13: 

d(dq)/dQ  *  «•••••*•••••••  V(q)  * _ *_  ft/s 

Form  14: 


Initial  Conditions  at  t®  »  days* 

Form  15: 

l_ _ l_ _ I _ !_• _ L* _ I _ L* _ _!_• _ I _ I _ I _ I _ I 

Form  18: 

Initial  Orbital  Central  Angle 
Form  20: 

dq(M)  *  . 

Form  21: 

dq(S)  ®  •••••••••••• 

Form  22: 

Initial  Lunar  Longitude 


Fora  23: 

Initial  Solar  Longitude 
Fora  24: 

Initial  Lunar  Perigee  Longitude 


Fora  30: 
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Orbital  Radius 

300000  n  mi 

Orbital  Inclination  to  tha  Ecliptic 

90,00  deg. 

Mean  Orbital  Rate 

•13142  rad/aolar  day 

Orbital  Period 

47.811  days 

Normal  Orbital  Acceleration 

,004220  ft/sec/sec 

Computing  Step  Size 

•500  days 

Initial  Orbital  Central  Angle 

20,0  deg 

Initial  Lunar  Longitude 

225,0  deg 

Initial  Solar  Longitude 

66,0  deg 

Initial  Lunar  Periges  Longitude 

180,0  deg 

Initial  Conditions  at  t=  ,0  days. 

dr/R  *  0 

dp  s 

0  n  mi 

d(dr/R)/dQ  *  0 

V(r)  ■ 

♦0  ft/s 

dq  (rad)  *  0 

dq  * 

,00  deg 

d(dq)/dQ  *  0 

V(q)  *  2774,4  ft/s 

dq(M)  *  0 

dq(S)  *  0 


II 
2| 
3| 
4| 
5| 
6| 
71 
3| 
9! 
10 1 
HI 
121 
13| 
14 1 
151 
161 
17| 

I 


3 

-25 

-134 

-366 

-750 

•1296 

-1996 

-2822 

-3733 

-4672 

-5577 

-6383 

-7032 

-7472 

-7667 

-7599 

-7273 


WFT 

ft/8 

0 

0 

-4 

-12 

-21 

—33 

-44 

-5»: 

-62 

-66 

-66 

-61 

-52 

-39 

-23 

-5 

14 

32 


a(rV 


.064 i 
‘•0059 
-.0158 
••0240 
.,0295 
••0316 
-.0300 
■•0248 
-.0164 
■•0056 
•0065 
*0190 
,0307 
•0405 
•0475 
•0509 
•0504 
•0458 


FTr)7«TvTqTri(qTTF(q'^  1  k"  TqThITqIS}  FouiTl 

1  ft/s  1  1  1 

1  1  1  1 

deg) 

1 

deg) 

1 

degl 

1 

deg  | 

1 

•00431 2774 l-.OOSSI -,00881  0|  6Sy  i32| 

-,0035 |2771| -.0087 |-,0089| 

ol 

239| 

67| 

117| 

- ,0114 |2768| - -0041 | - .0072 1 

01 

253) 

68| 

104) 

-.0185|2767|  .0042|-.0041| 

Ol 

266) 

69) 

93) 

-,0241|2769|  ,0151|-,0003| 

01 

280) 

70| 

84 1 

-.0279)27731  .02721  .00351 

01 

292) 

711 

78 1 

-.0295)27801  ,0393 |  .0072 | 

ol 

3051 

72| 

76| 

-.0290127901  .05061  .0108) 

ol 

317) 

73 1 

77 1 

-.026912802)  ,0602 |  .0144) 

ol 

329) 

74 1 

821 

-,0236|:816|  .0667)  .0174) 

ol 

341) 

75) 

881 

-.0199128321  .06911  .0193) 

11 

3531 

76) 

96| 

-.016212848)  .0667|  .0199) 

11 

5) 

77 1 

104| 

-.012912862)  ,05931  .0190 | 

11 

17) 

77| 

1111 

-.0106)20751  ,04721  .01681 

21 

29 1 

78| 

116| 

-.009612885)  .0315 |  .0136) 

2| 

41) 

79 1 

117| 

-,0102)28911  ,0134|  .0098) 

3| 

53) 

801 

1151 

-.0127) 2893I-.0052I  ,0059) 

3) 

66 1 

81) 

1091 

-.01701 2891 |-,0225|  .0024) 
1)11 

4| 

1 

78) 

1 

82) 

1 

100| 

1 

Initial,  Conditions  at  t*  17.0  days. 


dr/R 

* 

-2.42423402-02 

dr 

It 

-7273 

n  si 

d(dr/R)/dQ 

9 

1.14110033-02 

V(r) 

9 

31.7 

ft/s 

dq  (rad) 

m 

6,56216631-02 

dq 

9 

3.76 

<*•« 

d(dq)/dQ 

9 

6,80574022-02 

V(q) 

9 

2891.4 

ft/s 

dq(M)  *  -1,84920476-01 

dq(S)  «  -8,47662480-03 
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Orbital  Radius 

30000C  n  mi 

Orbital  Inclination  to  th«  Ecliptic 

90*00  dag* 

Naan  Orbital  Rata 

*13142  rad/solar  day 

Orbital  Pariod 

47*811  days 

Normal  Orbital  Accalaration 

*004220  ft/sac/sac 

Computing  Stap  Sine 

*500  days 

Initial  Orbital  Central  Angle 

20,0  dag 

Initial  Lunar  Longitude 

225.0  dag 

Initial  Solar  Longitude 

66*0  deg 

Initial  Lunar  Terigee  Longitude 

180.0  deg 

t  | 

wciomrsninTn  nmnnii  fiothi  m  rtomdi 

1 

days  | 

1 

n  ad  | ft/s j  |  Ift/sl  1  1 

!  1  i  1  1  1  1 

dag! 

1 

deg! 

! 

deg! 

1 

deg  | 

i 

“isl 

_r 

Kill 

83! 

89 1 

19 1 

-59551 

59 |  *0265|-*0309| 2879 |-*0467|-*0010| 

5! 

105! 

84! 

75J 

20 1 

-5063| 

66 |  *0132|-.0403i2870|-.0514i.,000s! 

51 

119! 

85 1 

611 

211 

-4104| 

68 !-*0022| -*0518) 2860 |-*0518|  *0003| 

6| 

133! 

86| 

47| 

22 1 

-3157| 

64| -*0210 | -*0668 | 2851 | -*0506|- *0021| 

6! 

147! 

87 1 

35) 

23| 

-2320| 

53i-*0409|-*0817| 2841 | -*0508 | -*0114 | 

6| 

1621 

88| 

261 

24| 

-1687| 

36 | -. 0454 | -. 0801 |2p?9|-. 0414 (-. 0145 | 

6| 

176! 

89| 

27 1 

25 1 

-12821 

21 | -*0347 | -.0653 | 2822 | -*0199 | -.0041 | 

6| 

191! 

90 1 

361 

26 1 

-1057| 

11 | -*0220 | -.0520 | 2820 ! -*0027 |  *0054| 

7| 

206! 

91| 

48| 

27 1 

-947| 

5 j -*01031 -.0420 | 2822 |  *0059|  *0097| 

7| 

220| 

92! 

611 

26| 

-891| 

3|  .00001 -.0342! 2825!  .0072 |  .0097 | 

7| 

234( 

93 1 

72! 

29 1 

-836| 

5|  ,0080| -.0283 |2828|  ,0034|  .0070! 

7| 

248! 

94 1 

81| 

30 1 

-740| 

9|  .0128|-,0245 | 2829 | -.0038 |  .00271 

7! 

262| 

95! 

88| 

311 

-580| 

14 |  .014ll-.0227l2827l-.0122t-.002l! 

7| 

275! 

96| 

91! 

32 1 

-349 1 

19 |  .0119 | -.02281 2824| -.0202 ! -.0065 | 

8! 

288] 

97! 

90 1 

33| 

-57| 

22|  .0068|-.0246i2818i-.0261t-,0099| 

8! 

301| 

98| 

87! 

34 1 

270| 

23 | -,0005| -.0278| 2811| -.0288) -.0118 I 

8! 

314 1 

98| 

81 1 

35 1 

593| 

22 I-.0089 I -.0318(2803! -.0275 |-.0il8| 

81 

326| 

99! 

72 1 

36| 

871| 

17 | -.0175 | -.0366 | 2797 | -.0216 !- .0095 | 

8! 

338| 

1001 

63! 

37| 

1058| 

9 |-.0254| -.0419 | 2792 | -.0113 | -.00481 

8| 

350! 

101! 

54 1 

38i 

1U4| 

-li-.0317|-.0476|279li  ,0027{  .0017! 

8| 

2! 

102! 

46! 

39 1 

1007| 

-14 | -.0350| -.0525|2794|  .0169 |  .0070! 

8! 

13! 

103! 

40) 

40| 

722! 

-26| -.03301 -.05 35 | 2799 (  .0258!  ,0068| 

S! 

25! 

104! 

39| 

41| 

267| 

-37l-.0256l-.0490! 28051  .0296!  .0026| 

81 

371 

1051 

43! 

42 1 

-319| 

-45|-.0161| -.0418128111  .032?!  .0000! 

8| 

49! 

106! 

5ll 

43| 

-989| 

-49l-.0062l-.0343! 2817!  .0357|-.0003| 

8| 

62! 

107! 

62! 

44 1 

L 

-1690! 

-49!  .0038!-.0268(2824!  ,0366|  .0000) 

III!!1 

9| 

! 

74 1 

_ L 

108| 

! 

74! 

1 

Initial  Conditions  at  t«  44,0  days. 


dr/R 

* 

-5.63352144-03 

dr 

■  -1690  n  mi 

d(dr/R)/dQ 

c 

-1,77619235-02 

Vir) 

*  -49.3  ft/e 

dq  (red) 

m 

1.48446267-01 

■  8.51  deg 

d(dq)/dQ 

9 

2.36641115-02 

V(q) 

*  2824.0  ft/s 

dq(N) 

9 

-1.81743453-01 

dq(S) 

9 

-2.34823162-02 
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Orbital  Radius 

Orbital  Inclination  to  tha  Ecliptic 

Ha  an  Orbital  Rata 

Orbital  Psriod 

Nornal  Orbital  Accalaration 

Computing  Stap  Sixa 

Initial  Orbital  Cantral  Angla 

Initial  Lunar  Longituda 

Initial  Solai*  Longituda 

Initial  Lunar  Parigaa  Longituda 


300000  n  ai 
90*00  dag* 

*13142  rad/solar  day 
47*811  days 
*004220  ft/sac/sac 
*500  days 
20*0  dag 
225*0  dag 
66,0  dag 
180,0  dag 


I  days  |  n  ai  lft/al  I  |ft/a 

1 _ I _ I _ I _ I _ I _ I 

|  5ST|  -2372 j  -46  i  Voi3^l-.0l9S|2i3Of 

I  46 |  -29831  -39 |  .0225 |-.0126| 2836! 

|  471  -34791  -30)  .02921 -.0066 | 2839 | 

i  48|  -3825|  -19 |  «0328|-*002l|2839i 

|  49 |  -40021  -7 |  *03221  .0004|2836|« 

|  S0|  -4015|  4|  *02721  *&004|2829|- 

II  I  I  I  I  I 

Initial  Conditions  at  t*  50*0  days* 


I  I 

_ I _ I 

*0  3&9 1  -*0b0&  I' 
*0271i-*0025| 
*0165|-*0061| 
*0031j-*0109| 
*0114 1 -.0163) 
.0248) -.0215! 

I  I 


dagj  dagj 

~rrw 

9|  100! 
9|  114| 
9|  128) 
10 1  1421 
10|  1561 
I  I 


dagj  dag  | 

109|  571 

110|  1001 
nil  1121 
112|  1221 
113)  1281 

1141  1291 

I _ I 


dr/R 

d(dr/R)/dQ 
tq.  (rad) 
d(dq)/dQ 
d*<IO 
d%(S) 


-1*33824949-02 

1.60698345-03 

1.74472449-01 

3*366235^7-02 

-1*26559426-01 

-2*68203919-02 


dr  *  -4915  n  ai 
V(r)  e  4.5  ft/a 
dg  *  10*00  dag 
t(<:>  *  2629.4  ft/s 


; 

i 

; 
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IDENTIFICATION  OF  JOSS  TERMS 

The  following  list  defines  the  terns  Appearing  in  the  preceding 

* 

; 

prograa  in  teros  of  the  variables  used  in  the  body  of  the  report 

wherever  possible. 

\ 

JOSS 

Definition 

i 

f 

r 

q 

&0  (red) 

t 

t  (days) 

i 

A(0) 

*0  (n  *i) 

i 

A(l) 

starting  value  o£  t 

f 

A<2) 

terminal  value  of  t 

r 

A(3> 

®0<de*> 

A(^) 

computing  step  size  (days) 

/ 

(J 

* 

A(5) 

printout  step  size  (days) 

; 

A(6) 

Sa(0) 

: 

A<7) 

8(0)  (deg) 

| 

A(8) 

A(9) 

n 

¥■ 

A(IO) 

ti/180  (rad  /day) 

r 

-V 

* 

A(U) 

aa(n  1) 

A(12> 

Vn  *l) 

: 

A(I3) 

c» 

[■ 

: 

A(14) 

E 

A(15) 

6^  (rad/deg) 

f 

A(16) 

ep  ("*> 

f 

A(17> 

rc  <ft) 

j  ; 

A(18) 

(tad/day) 

t  A' 

A(W) 

§  (rad /sec) 
c 

? 

A(20) 

unperturbed  orbital  velocity  (ft/sec) 

i  s 

A(21) 

unperturbed  orbital  acceleration  (ft/ sec") 

t 

A<22) 

unperturbed  orbital  period  (days) 

f  ‘ 

A/ 23) 

§  (rad /day) 

\  • 

A(24) 

80  (rad /day) 

i  ! 

A(25) 

'  *• 

? 

*  * 

A<26) 

$•5 

JOSS 

Definition 

A(29) 

printout  control  parameter 

A(30) 

termination  control  parameter 

A(31) 

^mo^o 

A(32) 

<Mo 

A(33) 

®rap/8o 

A(34) 

(Qjno  “  ®mp)/®o 

A(52) 

6mp(0)  (deg) 

A(53) 

®mp(0)  («d) 

A(61) 

0(0)  (deg) 

A(62) 

0(0)  (rad) 

i(l) 

6r  (n  mi) 

i(2) 

50  (deg) 

1(3) 

9m  (deg) 

i<4> 

®  (deg) 

t(5) 

ccs'1  (i  •  T  ) 
on 

1(6) 

d(6r)/dt  (ft/a) 

1(7) 

d2(5r/ro)/d0g 

1(8) 

a8r  +  amr 

1(9) 

r§ 

Jd) 

d2(60)/d6j 

j(2) 

•a0  +  ad0 

«d,  0) 

6r/r0 

1) 

d(5r/ro)/d0o 

z(2,  0) 

60  (rad) 

z(2,  1) 

d(60)/d0o 

r.(3,  0) 

60m  (rad) 

r(4,  0) 

69  (rad) 

Function* 


D(q, 

z) 

right  side  of  Eq.  (55) 

F(q, 

z) 

right  side  of  Eq.  (49) 

Mq, 

z) 

right  side  of  Eq.  (48) 

M(q, 

z) 

0m  (rad) 
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JOSS  Definition 

R(q»  z)  r/p 

S(q,  z)  0  (rad) 

a  Qf0  (rad) 

b(q,  z)  right  side  of  Eq.  (54) 

d(q)  («m/p)3 

e(q)  (aE/R)3 

h  computing  step  (rad) 

p(q»  z)  (p/rm)3 

r(z)  r/re 

x(M,  q,  z)  (i  •  im) 

x(S,  q,  z)  (I  •  rx) 

y(M,  q,  z)  (j  *  Im) 

y(S,  q,  z)  (j  •  ?i) 

It  should  be  emphasized  that  the  abo/e  definitions  apply  only  to 
the  program  for  Runge-Kutta  Integration  of  Orbital  Equations. 

PROGRAM  DESCRIPTION 
Part  1 

The  input  data  required  in  the  execution  of  this  part  specify 
the  following  aspects  of  the  problem: 

o  Reference  orbit 
o  Starting  time,  t^ 

o  Initial  geometry  at  t  *  0  (injection) 

o  Starting  values  of  6r/r0,  60,  6r'/r0,  68',  50  ,  6© 

m 

Part  2 

In  this  part,  the  parameters  necessary  in  the  solution  are 
evaluated. 

Part  3 

This  part  controls  the  printout  of  the  material  from  Orbital 
Radius  down  to  the  first  line  of  output  data  (see  sample  printout) . 
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Part  4 

This  part  results  in  the  printout  of  the  values  of  the  dependent 
variables  6r/r0,  60,  6r'/r0,  60  ',  <$9m,  and  6©  at  the  current  time,  t, 
as  well  as  6r(n,  mi),  66 (deg),  Vp  (ft/s)  and  Va  (ft/s). 

Part  5 

This  part  controls  the  integration  procedure,  the  printout  of 
the  data  at  the  end  of  each  printout  interval,  and  finally,  advances 
the  independent  variables  q  and  t,  after  which  the  procedure  is  re¬ 
peated  until  the  terminal  time  is  reached. 

Part  6 

This  part  computes  the  output  quantities  and  prints  them  as  one 
line  in  the  output  table  in  the  following  order:  t  (days),  6r  (n  mi), 
Vr  (ft/s),  d2(6r/ro)/d0g,  amv  +  asr,  Va  (ft/s),  60  (deg),  6m  (deg), 

@  (deg),  9sm  (deg). 

Part  7 

This  part  specifies  the  dependent  variable  to  be  evaluated  in 
the  integration  as  follows: 

k  Variable 

1  6r/r0 

2  60 

3  60m 

4  60 

Part  8 

This  part  specifies  the  order  of  the  derivative  of  the  dependent 
variable  by  the  value  of  j. 

Part  9 

Each  step  of  this  part  is  executed  for  6r/r0,  69,  fir7 /rQ,  60', 
69m,  and  6©  before  advancing  to  the  next  step.  Step  9.08  completes 
the  Runge-Kutta  Integration  over  one  computing,  interval  for  the  six 
quantities  listed  in  the  previous  sentence. 
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Orbital  Injection  and  Correction  Program 


Main  Control  Program 

********************** 


1*001  ^Initiation  (Do  part  1*) 

1*02  Demand  A(0)  as  "Reference  Orbital  Radius  (n  ml)". 

1*03  Demand  A(8)  as  "Initial  Orbital  Central  Angle  (deg)", 

1*04  Demand  A(6)  as  "Initial  Lunar  Longitude  (dag)"« 

1,05  Demand  A(7)  as  "Initial  Solar  Longitude  (deg)"* 

1,06  Demand  t  as  "Current  Tine  (days)", 

1,07  Demand  Z(1,0)  as  "Current  dr/R", 

1,08  Demand  Z(l,l)  as  "Current  d(dr/R)/dQ", 

1,09  Demand  Z(2,0)  as  "Currant  dq  (rad)", 

1.10  Dwstand  Z(2,l)  as  "Current  d(dq)/dQ", 

1.11  r«w®nd  Z( 3,0)  as  "Current  dq(M)  (rad)", 

1«12  Demand  Z(4,0)  as  "Current  dq(S)  (rad)", 

1.13  Type  "The  computation  can  be  shortened  if  the  reference  radius", 

1.14  Type  "is  the  sane  as  in  the  previous  case,", 

1.15  Line, 

1.16  B*»and  B  as  "The  radius  is  the  same  as  before,  (true  or  false)", 
1*17  Recall  item  12  (initA)  if  not  B, 

1.18  Recall  item  13  (initB)  if  B, 

1.19  To  step  1.25  if  B. 

1.20  Do  part  2, 

1*21  Do  part  3, 

1.22  Do  part  4  for  ns0(l)100, 

1.23  Do  part  16  for  i*2(l)7, 

1.24  Recall  item  13  (initB)  if  not  B, 

1.25  Do  part  8, 

1.26  Page, 

1.27  Do  part  13, 

1.28  Do  part  16  for  i«8(l)ll,13(l)15. 

1.29  Delete  all  forms, all  formulas, 

16,01  Delete  part  i. 


Item  12  (initA) 

Computation  of  Coefficients  in  Series  Expansion 

WmWmWWWWWWmnWmWWWWWnWmmVWMWWWWmWWWmWWIiWVWWWIlWivlV 


2,001  ^Evaluation  of  Parameters 
2,01  Set  m*82«27, 

2,02  Set  o*arg(-l,0)/180, 

2,03  Set  A(l)«238857*5280/6080, 

2,04  Set  A(2)»l/sqrt[(»-l)/n»A(l)*3/A(0)*3], 
2,05  Set  A(4)«,22998/A(2), 

2,06  Set  A(5)«.0172/A(4). 
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2.07  Set  K»2»A(0)»A(1)/[A(0)*24A(1)*23» 

2.08  Sat  V«A(4)»A(0)*6080/86400« 

2.09  Sat  k*[A(l)*K/2/A(0)3*1.5, 

2.10  Sat  C<0,0)«0. 

2.11  Delate  C. 

2.12  Sat  C(0,0)-«0. 

2.13  Let  C  be  sparse. 

2*14  Sat  b(0)*0, 

2.15  Let  b  be  sparse. 

3.001  ^Determination  of  P(i,j)  and  Q(i,j) 

3.01  Set  P(0,0)*-A<5)*2/4. 

3.02  Set  P(2,0)*-3«P<0,0). 

3.03  Set  P(0t2)*P<2,0). 

3.04  Set  P(2,2)=P<2,0)/2. 

3.05  Set  P(2,-2)=P(2,2). 

3.06  Set  Q(2,0)«-P(2,0). 

3.07  Set  Q(2,2)*-P<2,2). 

3.08  Set  Q(2,-2)«Q(2,2). 

3.09  Let  P  be  sparse. 

3.10  Let  Q  be  sparse. 

4.001  *Determination  of  b(n,i) 

4.002  Type  n  if  fp(n/10)*0. 

4.01  Set  b(0)»[r.*0:l}b(G)*(n-l)/n3  if  fp(n/2)*0. 

4,02  Set  b(l)=[nal:l5b(l)*n/(n+l)3  if  ffc(n/2)*0, 

4.03  Do  part  5  if  n*0  or  n*l. 

4.04  Done  if  n*0  or  n*l, 

4.05  Set  i«[fp(n/2)«0:2;33. 

4.06  Set  b(i)>«Cisn:l/2*(n-l);b(i-2)»(n-i+2)/(n+i)*[i*2t2;l33. 
4.07  Set  A(9)*diaj[i*20,i«19,i*n3. 

4.08  Do  part  5  if  A(9). 

4.09  Done  if  A(9). 

4.10  Set  isi+2, 

4.11  To  step  4.06. 

5,001  *Detenaination  of  C(i,j)  (Parts  5,6,7) 

5.01  Set  A(10)s[n*0il;A(10M2»n+l)»K/2/n3. 

5.02  Set  A(ll)*[fp(n/2)»0:0;13. 

5.03  Set  A(i2)*[fp(n/2)«0:20}193. 

5.04  Do  part  6  for  isA(ll)(2)A(12). 

6.01  Do  part  7  for  j*A(ll)(2)A(12). 

7.01  Set  C(i,j)*C(i,j)tk»A( 10)»b(i)»b(j)»[j*0:l}.5], 


p 
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Item  13  (inltB) 

Computation  of  Orbital  Injection  or  Correction 

********************************************** ft* 


8*001  ^Initialization  of  Simmations 

8*01  Set  A(16)*fpi:(A(6)t.22998^t/o)/3603»360t 

8.02  Set  A(17/«fp[(A(7)+t>Ci?2«t/o)/360]»360, 

8.03  Set  A(13)*fp[(A(8)+A(4)<>t/o+Z(2,0)/o)/3603*360, 

8.04  Set  A(21)*0c 
8.05  Set  A(22)s0. 

8.06  Set  A(23)*0. 

8.0?  Do  part  9  for  i«0(l)21. 

9.001  *Range  of  j  in  Sunationa 
9.01  Type  i  if  fp(i/5 )*0. 

9,02  Set  A(13)«[fp(i/2)a0:20;2l3. 

9,03  Do  part  10  for  j«-A(13)(2)A(13). 

10.001  *Lunar  Terms  in  Summations 
10.01  Set  A(21)«A(21)+F(p,q,v,u), 

10.02  Set  A(22)*A(22)-K3<p,q,v,u), 

10.03  Set  A(23)SA(23)+H(p,q,w(s)« 

10,04  Do  part  11  if  disj[conj(is0,0£ j£2),conj(i«2,-2£jS2)]. 
10.05  Done  if  di«j[i*21,j*2l3. 

10.051  *Deternination  of  Injection  Conditions 
10.06  Set  z(l,0)s2«A(21)-A<22). 

10.07  Set  z(l,l)«A(23). 

10,08  Set  z(2,0)*Z(2,0)« 

10.09  Set  z(2,l)e-3*A(21)+2*A(22). 

10.10  Set  z(3,0)«Z(3,0). 

10.11  Set  z(4,0)«Z(4,0). 

10*12  Set  A(24)*-Cp(0,0)+P(0f0)3«A(0)/3. 

10.13  Set  A(29)«z(2,l)-2»(Z<l,0)-z(l,0)). 

10.14  Set  A(30)*(1+Z(1,C))*(1+Z(2,1))»V, 

10.15  Set  A(31)a(l+Z(l,0))*(l+A(29))»V. 

10.16  Set  A(32)a2*A(0)-C2»(Z(l,0)-z(lt0))+Z(2,l)-z(2,l)3+A(24). 

10.17  Set  A(33)a-3*OZ(1,0)+Z(2,l)-A(21)3*A<4). 

10.18  Set  A(34)*C3»Z(1,0)+2*Z(2,1)-A(22)3*A<0). 

10.18  Set  A( 35)*[Z(1,1)-A(23)3*A(0). 

10.20  Set  A(36)*A(31)-A<30). 

10.21  Set  A(37)*(z(l,l)-Z(l,l))*V. 

10.22  Set  A( 38)*sqrt [ A( 36 )*2+A( 37 )*23 • 

10.23  Set  A( 39 )*sqrt[A( 34 )*2+A( 35 )*23 • 

10.24  Set  A(40)»2*A(39)/A(0)*A(4). 

10.25  Set  A(41)*(ltz(l>0))»<l+z(2,l))*V. 

10.26  Set  A(42)*z(l,l)*V. 

13.27  Sot  A(43)«sqrtCA(41)*2tz(l,l)*2»V*2l. 

10.28  Set  A(44)*arg[A(4l),z(lel)»V3/o. 

11,001  *Solar  Terms  in  Summations 
11,01  Set  A(21)sA(21)+F(P,Q,W,U). 

11,02  Set  A(22)«A(22)+G(P,Q,W,U), 
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11.03  Set  A(23)*A(23)+H(P,Q ,W,S), 

13.001  ^Printout 

13,01  Type  "Orbital  Inj*«tion"  if  t«0. 

13,02  Type  "Orbital  Correction"  if  t*0, 

13.03  Type  _,fora  1,_. 

13.04  Type  A(0)  in  fora  2. 

13.05  Type  fora  27, 

13.06  Type  V  in  fora  3. 

13,07  Type  A(6)  in  fona  5, 

13.08  Type  A(7)  in  fora  6. 

13.09  Type  A<8)  in  fora  4. 

13.10  Line. 

13.11  Typo  t  in  fora  34  +tv(t«»0). 

13.12  Do  part  14  if  t*0, 

13.13  Do  part  15  if  t*0. 

14.001  *Prin;out  of  Injection  Conditions 
14.01  Type  .fora  7, fora  24.  , 

14,02  Type  2i(l,0)«A(0),s(l,0T»A(0)  in  fora  8, 

14.03  Type  A<30),A(41)  in  fora  9, 

14.04  Type  Z(l,l)»V,z(l,l)«V  in  forn  10. 

14.05  Type  A(32),A(24)  in  fora  11. 

14.06  Type  A(33),0  in  fora  12. 

14.07  Type  A(39),0  in  fora  13. 

14.08  Type  A(40).0  in  fora  14. 

14.09  Line. 

14.10  Type  A(43)  in  fora  30. 

14.11  Type  A(44)  in  fora  31. 

14.12  Type  ,  .fora  25tfora  26,  . 

14.13  Type  7(T.0),z(l,0)  in  fora  18. 

14.14  Type  Z(l,l),z(l,l)  in  fora  19. 

14.15  Type  Z(2.0),z(2.0)  in  fora  20. 

14.16  Type  Z(2,l),z(2,l)  in  fora  21. 

14.17  Type  Z(3,0),z(3,0)  in  fora  22. 

14.18  Type  Z(4,0),s(4,0)  in  form  23. 

15.001  ^Printout  of  Correction  Conditions 
15.01  Type  .fora  7, fora  24,  . 

15.02  Type  ?(1,0)*A(0),Z(1,0T*A(0)  in  fora  8. 

15.03  Type  A(30).A(31)  in  fora  9. 

15.04  Type  Z(l,i).V,*(l,l)*V  in  fora  10e 
15.05  Type  A(32),A(24)  in  fora  11. 

15.06  Type  A(33),0  in  fora  12. 

15.07  Type  A(39),(Z(1.0)-z(l,0))»A(0)  in  fcra  13. 
15.08  Type  A(40),2*(Z(l,0)~z(l,0))*A(4)  in  fora  14. 
15.09  Line. 

15.10  Type  A(36)  in  fora  15. 

15*11  Type  A(37)  in  fora  16, 

15.12  Type  «(36)  in  fora  17. 

15.13  Type  ,  .fora  25, fora  26,  , 

15.14  Type  7(T,0),Z(1,0)  in  fora  18, 

15.15  Type  Z(l,l),z(l,l)  in  fora  19. 

15.16  Type  Z(2,0),Z(2,0)  in  fora  20, 


15.17  Type  Z(2,1),A(29)  in  fo»  21, 

15.18  Type  Z(3,0),Z(3,0)  in  form  22, 

15.19  Type  Z(4,0),Z(4,0>  in  form  23, 


F(p»<l»w»u):  -l/3»[|i|+| j |»0:2*p(0,0);3»q(i ,j)»u(i,j)/w(i,j)3 
(p(i,j)-2*q(i,j)*w(i, j))/(w(i,j)*2-l)»u(i,j ) 

H(p»q»w»u):  (p(i,j)*w(itj)-2»q(i,j)>u(i,j)/(w(i,j)A2~l) 

S(i,j):  sin(h(i,j)) 

U(i,j):  cos(h(i,j)) 

W(i,j):  l+j«A(5) 

c(i,j)i  C(i,|j|) 

d(i,j):  -A(i)/A(0)/4*[sura[r*|  j|+i,|  j|-l:c(i-l,r)-c(i+l,r)3-l(i»j)3 

A(l)/A(0)/4*Csun[rBi+l,i-l:s«t«C8s|j  1+1,1  j|-l:c(r,s)]]-l(i,j)l 
g(i,j):  fp((i»A(18)+j«A(16))/360)*360»o 

h(i,j):  fp((i*A(18)+j  *A(17))/360)«360»o 

l(i, j):  C|i’j!*l:2;0]-[i«ljc<0,|j|+l)+c(0,|j|-l)j0] 

p(i,j):  C-c(i,j)+e(i,j)].A(2)*2/» 

q(i, j):  d(i,j).A(2)*2/n 

sin(g(i, j )) 
u(i,j):  coe(g(i,j)) 
w(i,j):  i+j »A(2) 


For*  1: 

Reference  Condition 
Fora  2: 

Reference  Orbital  Radius  (n  mi) 


Fora  3: 

Reference  Orbital  Velocity  (ft/s)  , 

Fora  4: 

Initial  Orbital  Central  Angle  (deg)  _ , 

Fora  5: 

Initial  Lunar  Longitude  (deg)  _ , _ 

Fora  6: 

Initial  Solar  Longitude  (deg)  _ • 

Fora  7; 

Current  Corrected 


Fora  8: 

Increment  in  Orbital  Radius  (n  mi) 
Fora  9: 

Horizontal  Velocity  (ft/s) 


Fora  10 1 

Radial  Velocity  (ft/s) 


Fora  11: 

Average  Radial  Increment  (n  mi) 
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Fora  12: 

Average  Orbital  Rata  Incramant  (rad/ day) 

Form  13: 

Oscillatory  Amplitude  in  dr  (n  mi)  _ 

Form  14: 

Oscillatory  Amplitude  in  d(dq)/dt  (rad /day) 

Form  IS: 

Horizontal  Velocity  Correction  (ft/s) 

Form  16: 

Radial  Velocity  Correction  (ft/s) 

For*  17: 

Total  Velocity  Correction  (ft/s) 

Form  18: 

dr/R  . . .  . . 

Form  19: 

d(dr/R)/dQ  .  . . . . 

Form  20: 

<*q  . . 

Form  21: 

d(dq)/dQ  . . . . 

Form  22: 

dq(M)  . . . . 

Form  23: 

dq(S)  . . . 

Form  24: 

Form  25: 

Currant  Corrected 

Form  26: 

Value  Value 

Form  27; 

Orbit il  Inclination  to  Ecliptic  (deg) 


Value 


Value 


90 


Form  30: 

Injection  Velocity  (ft/s) 
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For*  31: 

Injection  Peth  Angle  (deg) 

Font  34: 

Injection  et  _  day* 

Fom  35: 

Correction  et  •  .  days 
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Orbital  Injection 
Reference  Condition 


Reference  Orbital  Radius  (n  mi)  303444 
Orbital  Inclination  to  Ecliptic  (deg)  90 
Reference  Orbital  Velocity  (ft/s)  2758*59 
Initial  lunar  Longitude  (deg)  158*80 
Initial  Solar  Longitude  (deg)  *00 
Initial  Orbital  Central  Angle  (deg)  *00 


Inje*  tion  at  0  days 


Current  Corrected 

Value  Value 

Xr.creaient  in  Orbital  Radius  (n  ml)  *00  2318*40 

Horizontal  Velocity  (ft/e)  2758*59  2790*09 

Redial  Velocity  (ft/e)  *00  9*80 

Average  Radial  Increaent  (r<  mi)  ->10021*88  1528*93 

Average  Orbital  Rate  Increment  (rad /day)  *0073751  *0000000 

Oscillatory  Amplitude  in  dr  (n  ml)  9293*12  *00 

Oscillatory  Amplitude  in  d(dq)/dt  (rad /day)  *0079128  *0000000 

Injection  Velocity  (ft/s)  2790*11 

Injection  Path  Angle  (deg)  *201 


Current  Corrected 

Value  Value 

dr/ R  0  7*64029580-03 

d(dr/R)/dQ  0  3.55243699-03 

dq  0  -7 

d(dq)/dQ  0  3*74893350-03 

dq(H)  0  0 

dq(S)  0  0 


Orbital  Correction 
Ref* rence  u»di.tio» 
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Reference  Orbital  Radius  (n  ai) 

303888 

CTbitel  Inclination  to  Ecliptic  (deg) 

90 

Reference  Orbital  Velocity  (ft/s) 

2758.59 

Initial  Lunar  Longitude  (deg) 

158.80 

Initial  Solar  Longitude  (deg) 

•CO 

Initial  Orbital  Central  Angle  (deg) 

.00 

Correction  at  235.5  days 

Current 

Corrected 

Value 

Value 

In cv saint  in  Orbital  Radius  (n  ai) 

1986.31 

1986.31 

Horizontal  velocity  (ft/s) 

2771.38 

2603.57 

Radial  Velocity  (ft/s) 

-18.78 

-75.48 

Averse  Radial  Inersseant  (n  ai) 

-5516.00 

1526*93 

Average  Orbital  lata  Increasnt  (rad/day) 

.0088976 

*0000000 

Oscillatory  Aaplitude  in  dr  (n  ai) 

9925*58 

20*28 

Oscillatory  Aaplituds  is  d(dq)/dt  (red /day) 

.0080256 

*0000172 

Horisoctel  Velocity  Correction  (ft/s) 

32*22 

Radial  Velocity  Correction  (ft/a) 

-56.78 

Total  Velocity  Correction  (ft/s) 

65,25 

Current  Corrected 

Talus  Value 


<WR 

d(dr/R)/dQ 

d(dq)/<5Q 

dqOO 

4q($) 


{.M59000M3 

-€.7*925278-03 

6*22908691-01 

-1.9.1080*92-03 

7.18890938-02 

-5.9*390258-02 


6*58590009-03 

-2.73623055-02 

6*22808691-01 

9*69817132-03 

7*18880936-02 

-5*92090258-02 
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IDENTIFICATION  OF  JOSS  TERMS 

The  following  list  defines  the  terms  appearing  in  the  preceding 
program  in  terms  of  the  variables  used  in  the  body  of  the  report 
wherever  possible. 


JOSS 

Definition 

K 

K 

V 

vo 

/Kdm\3/  2 

k 

l  2rJ 

m 

U 

n 

n 

o 

tt/180  rad /deg 

t 

ti  (days) 

A(0) 

r0  (n  mi) 

A(l) 

am  (n  mi) 

A(2) 

^mo/®o 

A(4) 

§0  (rad /day) 

A(5) 

®o^®o 

A(6) 

0m(O)  (deg) 

A<7) 

0(0)  (deg) 

A(8) 

0(0)  (deg) 

A(16) 

0m(ti)  (deg) 

A(17) 

9(ti)  (deg) 

A(18) 

9(t!>  (deg) 

A(21) 

B0/3 

A(22) 

Al 

A(23) 

-a2 

A(24) 

6r8E  (corrected)  (n  mi) 

A(29) 

66' 

A(30) 

Vh  (current)  (ft/ sec) 

A(31) 

Vh  (corrected)  (ft/sec) 

A(32) 

6rss  (current)  (n  mi) 

A(33) 

60ss  (current)  (rud/day) 

A(36) 

AVjj  (ft/sec) 

A<37) 

AVp  (ft/sec) 

v 
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JOSS 

Definition 

A(38) 

AV  (ft/sec) 

A(39) 

6rogc  (current)  (n  mi) 

A(40) 

60Q8C  (current)  (rad/day) 

A(41) 

VH  (injection)  (ft/sec) 

A(42) 

Vr  (injection)  (ft/sec) 

A(43) 

V  (injection)  (ft/sec) 

A(44) 

Yp  (injection)  (deg) 

b(i) 

bni 

C(i»  j) 

Cij 

P(i,  J) 

Pij 

Q(i,  j) 

Qij 

Z(l,  0) 

6ri/r0 

2(1,  1) 

6ri/r0 

Z<2,  0) 

60i  (rad) 

Z(2,  1) 

60{ 

Z(3,  0) 

60m  (rad) 

Z(4,  0) 

6©  (rad) 

*(1»  0) 

6rp/r0 

«<l,  1) 

6jTD/r0 

z(2,  0) 

60p  (rad) 

z(2,  1) 

68d 

2(3,  0) 

60m  (rad) 

2(4,  0) 

60  (rad) 

Functions 

F(P,q,w,u) 

ij  increment  to  B0/3 

G(p,q,w,u) 

ij  increment  to  Ai 

K(p,q,w,8) 

ij  increment  to  -A2 

S(i,j) 

sin  (i0(tx)  +  j©(t!)) 

U(i,j) 

cos  (i0 (tx)  +  j0(tL)) 

W(i,j) 

nij 

c(i,j) 

Cij 

g(i,j) 

i9(ti)  +  j^mCtl)  (rad) 

h(i,j) 

i9(t!)  +  j© (tx)  (rad) 

JOSS 

Definition 

P(i»j) 

pu 

q(i.j) 

-u 

s(i, j) 

sin 

<ie(t!) 

u(i,j) 

cos 

(ie<tl) 

+  j®m(tl)) 

w(i,j) 

As  before,  these  definitions  apply  only  to  the  Orbital  Injection 
and  Correction  Program. 

PROGRAM  DESCRIPTION 
Part  1 

The  input  data  required  in  the  execution  of  this  part  specify 
the  following  aspects  of  the  problem: 

o  Reference  orbit 

o  Initial  geometry  (t  =  0)  (injection) 
o  Current  time,  t^ 

o  Current  values  of  the  dependent  variables,  6r^/r0,  69^, 

fir.'/r  ,  60,',  (59  ,  and  6© 

1  o’  1’  m 

Next,  Item  12  (initA),  which  includes  parts  2  through  7,  is  executed 
if  the  values  of  C^j  corresponding  to  this  referen.  .  orbit  are  not 
already  available  in  storage;  otherwise  the  control  is  transferred 
directly  to  Item  13  (initB),  which  includes  parts  8  through  18. 

Item  12  (initA) 

This  subroutine  evaluates  the  coefficients  Cj_j  in  the  series 
expansion  of  (am/rQ)^. 

Part  2 

In  this  part,  the  parameters  necessary  in  the  solution  are 
evaluated. 

Part  3 

This  part  evaluated  the  coefficients  P^j  and  Q^j  in  the  expan¬ 
sion  of  the  solar  perturbing  functions  agr  and  a^g,  respectively. 


**{  p.-y^  fjr' 
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Part  4 

This  part  evaluates  the  coefficients  bn^  of  the  Fourier  expan¬ 
sion  of 


i 

max 


where  n  takes  on  all  integer  values  from  0  to  100,  while  i  is  equal 

°  max 

to  20  or  n,  whichever  is  smaller.  After  this  part  is  completed  for  a 
given  value  of  n,  parts  5,  6,  and  7  are  executed  before  advancing  to 
the  next  value  ol  n. 


par  t  5 

This  part  determines  the  values  of  i  considered  in  the  evalua¬ 
tion  of  C^j. 

Part  6 

This  part  determines  the  values  of  j  considered  in  the  evalua¬ 
tion  of  Cj[j. 

Parr  7 

This  part  evaluates  the  contribution  of  the  nth  term  of  the 
binomial  expansion  of  (sm/t0)^  to  eac  ;  of  Cj:  j  elements. 


Item 


This  subroutine  determines  the  orbitaL  injection  conditions  if 
t  *  0,  and  the  orbital  correction  conditions  if  t  ^  0. 


Part  8 

This  part  determines  the  current  geometry  at  time  t^  by  evalua¬ 
ting  QCt!),  SmCti)  and  8(t^)  and  initiates  the  summations  involved 
in  determining  A^,  A2  and  B0. 
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Part  9 

This  part  specifies  the  values  of  j  included  in  the  summations 
for  Aj,  A2,  B0. 

Part  10 

In  this  part,  the  values  of  A2,  and  B^  are  computed  and  then 
used  to  determine  the  quantities  dr^/r^  dr^/ro,  and  60^  after  which 
various  other  output  quantities  are  computed. 

Part  11 

This  part  computes  the  contributions  of  the  solar  terms  to  Ax, 

A2>  and  Bq. 

Part  13 

This  part  controls  the  printout  through  the  value  of  the  Initial 
Orbital  Central  Angle  (see  sample  printout) ,  after  which  control  is 
transferred  to  part  14  if  orbital  injection  is  being  determined 
(tL  -  0),  or  to  part  15  if  an  orbital  correction  is  required  +  0) . 

Part  14 

This  part  controls  the  printout  of  the  orbital  injection  condi- 
tions  as  follows: 

o  The  current  and  corrected  values  of  dr  (n  mi),  V_  (ft/sec), 
(ft/sec),  drg8  (n  mi),  d0##  (rad/day),  dro;  (n  mi),  and 
Q8C  (rad/day). 

o  The  orbital  injection  velocity,  V  (ft/sec),  and  the  orbital 
injection  path  angle,  (deg)* 

o  The  current  and  corrected  values  of  dr/ro,  dr'/r  ,  d6,  d97 , 
d0m,  and  SB  in  the  form  needed  as  input  to  the  Runge-Kutta 
integration  program. 

Part  15 

This  part  controls  the  printout  of  the  orbital  correction  condi¬ 
tions.  It  differs  from  part  14  in  printing  the  horizontal  and  radial 
components  of  the  correction  velocity,  AVg  and  AVg,  as  well  as  its 

total  magnitude  AV,  all  in  ft/sec  rather  than  V  and  y  . 

P 


references 
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